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The role of species interactions in controlling the interplay between the stability of an ecosystem 
and its biodiversity is still not well understood. The ability of ecological communities to recover 
after a small perturbation of the species abundances (local asymptotic stability) has been well 
studied, whereas the likelihood of a community to persist when the interactions are altered (struc¬ 
tural stability) has received much less attention. Our goal is to understand the effects of diversity, 
interaction strenghts and ecological network structure on the volume of parameter space leading 
to feasible equilibria, i.e., ones in which all populations have positive abundances. We develop a 
geometrical framework to study the range of conditions necessary for feasible coexistence in both 
mutualistic and consumer-resource systems. Using analytical and numerical methods, we show 
that feasibility is determined by just a handful of quantities describing the interactions, yielding 
a nontrivial complexity-feasibility relationship. Analyzing more than 100 empirical networks, we 
show that the range of coexistence conditions in mutualistic systems can be analytically predicted 
by means of a null model of random interactions, whereas food webs are characterized by smaller 
coexistence domains than those expected by chance. Finally, we characterize the geometric shape 
of the feasibility domain, thereby identifying the direction of perturbations that are more likely to 
cause extinctions. Interestingly, the structure of mutualistic interactions leads to very heterogeneous 
responses to perturbations, making those systems more fragile than expected by chance. 

Natural populations are faced with constantly varying environmental conditions. Environmental conditions affect 
physiological parameters (e.g., metabolic rates m) as well as ecological ones (e.g., the presence and strength of 
interactions between populations mg). Therefore in order to persist, ecological communities necessarily need, at the 
very least, to be able to cope with small environmental changes. Mathematically, this translates into an argument 
on the robustness of the qualitative behavior of an ecological dynamical system: to guarantee robust coexistence, a 
model describing an ecological community needs at least to be (qualitatively) insensitive to small perturbations of 
the parameters [3 [7]. This notion has been formalized in the measure of robustness [8] or “structural stability” [9], 
expressed as the volume of the parameter space resulting in the coexistence of all populations in a community. 

While the local asymptotic stability (the ability to recover after a small change in the population abundances) 
of ecological communities has been studied in small [10] and large HTHni systems, the study of structural stability 
(i.e., the ability of a community to retain the same dynamical behavior if conditions are slightly altered)—despite 
being proposed early on as a key feature in the context of the diversity-stability debate [TMl8] —has historically been 
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restricted to the case of small communities, with the first studies of larger communities appearing only recently isms], 
and—because of mathematical limitations—dealing exclusively with the case of large mutualistic communities. Studies 
of structural stability have so far focused on the effect of ecological network structure (who interacts with whom) on 
the volume of parameter space leading to feasible equilibria, in which all populations have positive abundances. 

Here we develop a geometrical framework for studying the feasibility of large ecological communities. We overcome 
the limitations that have hitherto prevented the study of consumer-resource networks, thereby providing a unified 
view of feasibility in ecological systems. Using a random matrix approach (which helped identify main drivers of local 
asymptotic stability), we pinpoint the key quantities controlling the volume of parameter space leading to feasible 
communities, as well as its sensitivity to changes in these parameters. We then contrast these expectations for 
randomly connected systems with simulations on structured empirical networks, quantifying the effects of network 
structure on feasibility. 


Theoretical framework 


For simplicity, we consider a community composed of S species whose dynamics is determined by a system of 
autonomous ordinary differential equations: 


drii 

dt 


= Ui 



s 

i=i 



( 1 ) 


where rii is the density of population i, ri is its intrinsic growth rate, and Aij (which in principle could depend on n; 
see Supplementary Information) measures the interaction strength between population i and j. A fixed point n* (i.e., 
a vector of densities making the right side of each equation zero) is feasible if n* > 0 for every population. A fixed 
point is locally asymptotically stable if, following any sufficiently small perturbation of the densities, the system returns 
to a small vicinity of the fixed point. The fixed point is globally asymptotically stable if the system eventually return 
to it, starting from any positive initial condition within a finite domain. A system with a fixed point is structurally 
stable if, following a sufficiently small change in the growth rates r^, the new fixed point is still feasible and stable. 

To study the range of conditions leading to stable coexistence, we need to disentangle feasibility and local stability. 
This problem is well discussed in Rohr et al. [5], where it was solved for the case of one possible parameterization 
of mutualistic interactions. If A is diagonally stable or Volterra dissipative (i.e., there exists a positive diagonal 
matrix D such that DA -f A^D is stable), then any feasible fixed point is globally stable [201 HI]- Unfortunately, 
a general characterization of this class of matrices is unknown |22j . We proceeded therefore by considering only the 
matrices such that all the eigenvalues of A -k A^ are negative (i.e., the matrix A is negative definite in a generalized 
sense [ 23 ] j corresponding to D being equal to the identity matrix; see Methods). This choice reduces the number of 
parameterizations one can analyze, as not all the diagonally stable matrices are negative definite. However, as shown 
in the Supplementary Information, only very few parameter combinations are excluded from this set. Moreover, the 
effects of negative definitness are well-studied for random matrices |24] , and by using it we can extend the study of 
feasibility to any ecological network, including food webs. 

Our goal is to measure the fraction of growth rate combinations, out of all possible combinations, that lead to the 
coexistence of all S populations. Since we can separate stability and feasibility, we only need to find those leading 
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to feasible fixed points, and the condition above ensures that these will be globally stable. As pointed out before [S], 
the problem is not to find a particular set of ri leading to coexistence, but rather to measure how flexibly one may 
choose these rates. As shown in Fig. this quantity-indicated by S henceforth-can be thought of as a volume, or 
more precisely a solid angle, in the space of growth rates (see Methods). 

To calculate S, one might naively wish to perform direct numerical computation of the fraction of growth rates 
leading to a feasible equilibrium. While a direct calculation is viable when S is sufficiently small, this procedure 
becomes extremely inefficient for large 5" [9]. In the Supplementary Information we introduce a method that can be 
used to efficiently calculate S with arbitrary precision, even for large S. Using this method, we can accurately measure 
the size of the feasibility domain, with larger values of S corresponding to larger proportions of conditions (intrinsic 
growth rates) compatible with stable coexistence. For reference, we normalize S so that S = 1 when populations are 
self-regulated and not interacting (Methods), i.e., when the interaction matrix A is a negative diagonal matrix, and 
thus Eq. [l] simplifies to S independent logistic equations. 


Results 

Feasibility is universal for random matrices. 

May’s seminal work m pioneered the use of random matrices as a reference, or null model, of ecological interactions. 
A particularly interesting feature of random matrices is that the distribution of their eigenvalues (determining local 
stability) is universal This means that local stability depends on just a few, coarse-grained properties of the 
matrix (i.e., the number of species and the first few moments of the distribution of interaction strengths) and not 
on the finer details (e.g., the particular distribution of interaction strengths; see Supplementary Information). In 
fact, these moments can be combined into just three parameters: Ei, E 2 , and Ec (Methods). Together with S, they 
completely determine local asymptotic stability. 

We tested whether universality also applies to feasibility. We considered different random matrix ensembles obtained 
for different connectance values and distributions from which the matrix entries were drawn, but with constant values 
of S and of Ei, E 2 , and Ec- We then checked whether the size S of the feasibility domain depended only on 
these four quantities or also on finer details. Surprisingly, we found that the feasibility of random matrices is also 
universal (Methods and Supplementary Information). Two very different (random) ecosystems, with completely 
different interaction types and distributions of interaction strengths, but having the same number of species S and the 
same Ei, E 2 , and Ac, have the same S in the large S limit. This result has important theoretical implications, as it 
indicates those moments as the drivers of feasibility, but also very practical consequences, namely that the parameter 
space one needs to explore is dramatically reduced. 


An analytical complexity—feasibility relationship 

The universality of S suggests that it is amenable to analytical treatment. As explained in the Supplementary 
Information and shown in Fig.[^ when the mean and variance of interaction strengths are not too large (Supplementary 
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Information), we are able to derive the following approximation for ^ for random interaction matrices A: 

where S is is the number of species, d is the mean of A’s diagonal entries, and Ei = C/r, the product of the connectance 
C and the average interaction strength /r (see Methods). A more accurate formula is presented in the Supplementary 
Information. 

In analogy with the celebrated result of May m connecting stability and complexity, Eq. can be considered as a 
complexity-feasibility relationship. While in May’s scenario and in its generalizations [12j the effect of complexity and 
diversity on stability is always detrimental, it does depend on the interaction type in the case of feasibility. Given that 
d is negative by construction (Supplementary Information), having more species or connections can either increase 
{El > 0) or shrink {Ei < 0) the size of the feasibility domain, as a function of the sign of interaction strenghts (see 
Fig.§. It is important to stress that we computed S under the assumption of A being negative definite. When we 
consider how S depends on S and other parameters, we need to take into account the conditions making the matrix 
negative definite (see Methods and Supplementary Information). In the case of positive interaction strengths, this 
condition is d + SCp < 0, implying an upper bound for p, that depends on S. 


Our analytical formula predicts feasibility of empirical mutualistic networks and overestimates that of food 

webs 

Having explored the feasibility of random networks, we proceed to investigate the effects of incorporating empirical 
network structure. Ecological networks are in fact non-random |27H29j . and many studies have hypothesized that 
the structure of interactions could increase the likelihood of coexistence [501132) . Having an analytical prediction for 
random matrices, we can study whether it predicts the size of the feasibility domain for empirical networks as well. 
Fig. I shows the simulated values of for 89 mutualistic networks and 15 food webs (Supplementary Information), 
parameterized multiple times and compared with our analytical approximation (see Methods). We find that S of 
empirical mutualistic networks is well predicted by our formula, while it overestimates the feasibility domain of food 
webs, indicating that their non-random structure has a strong negative effect on feasibility. 

In the Supplementary Information, we compare the effect of the empirical structure of mutualistic networks with 
randomizations, by controlling for the interaction strengths. We show that, in the absence of variability in interaction 
strengths, the structure of empirical mutualistic networks has a positive effect on feasibility, which is strongly reduced 
when interaction strengths are allowed to vary. While this effect of empirical mutualistic networks is statistically 
significant, its effect on S is negligible compared to the effect of mean interaction strengths, and can only be detected 
by controlling very precisely for interaction strengths (see Supplementary Information). On a broader scale, as shown 
in Fig. the size of the feasibility domain of empirical networks is well predicted by our analytical formula. 

On the other hand, the negative effect of food web structure on S is substantial. In the Supplementary Information 
we compare each network with randomizations and also with predictions of the cascade model m. which has recently 
been shown to predict well the stability of empirical food webs [14] . By analyzing different parameterizations we 
found that the feasibility domain of empirical structures is consistently and significantly smaller than that of both 
the randomizations and the cascade model. For most of the webs, the prediction obtained from the cascade model 
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is better than that of randomizations, suggesting that the directionality of empirical webs plays a role in reducing 
feasibility, with other properties of the structure of empirical networks also contributing significantly to feasibility. 


The shape of the feasibility domain carries information on the response to perturbations, and can be 

analytically predicted for random interactions 

So far, we have focused on the volume of the parameter space resulting in feasiblity. However, two systems having 
the same S can still have very different responses to parameter perturbations, just as two triangles having the same 
area need not to have sides of the same length (Fig. [^. The two extreme cases correspond to a) an isotropic system 
in which, if we start at the barycenter of the feasibility domain, moving in any direction yields roughly the same 
effect (equivalent to an equilateral triangle); b) anisotropic systems, in which the feasibility domain is much narrower 
in certain directions than in others (as in a scalene triangle). For our problem, the domain of growth rates leading 
to coexistence is—once the growth rates are normalized—the {S — l)-dimensional generalization of a triangle on a 
hypersphere. For S = 3, this domain is indeed a triangle lying on a sphere as shown in Fig.[^ If all the S{S—l)/2 sides 
of this (hyper-)triangle are about the same length, then different perturbations will have similar effects on the system. 
On the other hand, if some sides are much shorter than others, then there will be changes of conditions which will 
more likely impact coexistence than others. We therefore consider a measure of the heterogeneity in the distribution of 
the side lengths (Fig.[^and Supplementary Information). The larger the variance of this distribution, the more likely 
it is that certain perturbations can destroy coexistence, even when S is large and the perturbation small. This way of 
measuring heterogeneity is particularly convenient because it is independent of the initial conditions. Moreover, the 
length of each side can be directly related to the similarity between the corresponding pair of species (Supplementary 
Information), drawing a strong connection between the parameter space allowing for coexistence and the phenotypic 
space. As in the case of S, this measure is a function of the interaction matrix and corresponds to a geometrical 
property of the coexistence domain. 

While S is a universal quantity for random networks, the distribution of side lengths is not: it depends on the full 
distribution of interaction strengths (Supplementary Information). On the other hand, as shown in the Supplementary 
Information, it is possible to compute it analytically in full generality, i.e., for any distribution of interaction strengths 
and any interaction types. In particular we are able to obtain an expression for its mean and variance, which depend 
only on S, Ei, E 2 , and Ec (Supplementary Information). Fig. shows that the analytical formula, in the case of 
random A, matches the observed mean and variance of side lengths of random networks perfectly. 

Empirical network structures correspond to more heterogeneous shapes 

As we have done for S, we can now test how non-random empirical network topologies influence the distribution of 
side lengths. Fig. shows that empirical food webs and, in particular, empirical mutualistic networks display a much 
larger variation in side lengths than expected by chance. This result is particularly relevant, indicating that even if 
the feasibility domains of empirical mutualistic networks are larger than those of random networks, their shapes are 
less regular than expected by chance, and thus we expect perturbations in certain directions to quickly lead out of 
the feasible domain of growth rates. 
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Discussion 

A classic problem in mathematical ecology is determining the response of systems to perturbations of model pa¬ 
rameters. In the community context, one important application is getting at the range of parameters allowing for 
species coexistence Several methods exist for evaluating this range HHMH], but they either rely on raw 

numerical techniques, or else can only evaluate system response to small parameter perturbations. Here, in the 
context of the general Lotka-Volterra model, we have given a method for the global assessment of all combinations 
of species’ intrinsic growth rates compatible with coexistence—what we have called the domain of feasibility. Our 
geometrical approach can determine not only the total size of the feasibility domain, but also its shape: it is always a 
simply connected domain forming a convex polyhedral cone whose side lengths can be evaluated from the interaction 
matrix. Applying our method to empirical interaction networks, we were able to characterize the region of parameter 
space compatible with coexistence; the importance of this kind of information is underlined by a rapidly changing 
environment that is expected to cause substantial shifts in the parameters influencing these systems. 

The geometrical framework we employed, pioneered by Svirezhev and Logofet |25j . allows for the formulation of a 
complexity-feasibility relationship. In analogy with the celebrated complexity-stability relationship, it relates the size 
of the feasibility domain with diversity, connectance and interaction strengths of a random interacting community. 
While communities are not random, this relationship sets a null expectation for the scaling of the proportion of 
feasible conditions. We obtain that the mean of interaction strengths sets the behavior of feasibility with the number 
of species. If the mean is negative (e.g., in case of competition or predation with limited efhciency), the larger the 
system is, the smaller is the set of conditions leading to coexistence, while for positive mean (e.g., in the case of 
mutualism) the converse is true. 

Several recent works have studied the effect of network structure on coexistence in species-rich communities, with 
contrasting results [9l 1301 - 1321 [39] . Here we have shown that the fraction of conditions compatible with coexistence 
is mainly determined by the number and the mean strength of interactions. In terms of network properties, the 
relevant quantity is the connectance, with other properties (e.g., nestedness or degree distribution) having minimal 
effects. In particular, once the connectance and mean interaction strength are fixed, the matrices built using empirical 
mutualistic networks have feasibility domains very similar to that expected for the random case, as was also observed 
previously in a similar context |40) . 

The empirical network structure of mutualistic networks has a statistically significant effect on the size of the 
feasibility domain (see Supplementary Information). Whether this effect is ecologically relevant depends on the 
specific application at hand. For instance, the effect of structure could be neglected to quantify how the feasibility 
domain would change if a fraction of pollinators went extinct, and it could be evaluated using our analytical result. 
In contexts where the interaction strengths are strongly constrained, structure would play an important role. Our 
method provides, in this respect, a direct way of quantifying the importance of different factors, disentangling the 
way different interaction properties affect feasibility. 

For mutualistic interaction networks, our results clearly show which properties determine the global health of the 
community, and therefore indicate which properties should be measured in the field. While not observing a link or 
measuring a wrong interaction coefficient could have strong effects on ecosystem dynamics, they have very little effect 
on how the community copes with environmental perturbations and how likely extinctions are |41] . The major role 
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is played by corse-grained statistical properties of the interactions, such as connectance or the mean and variance of 
the interaction strengths. 

For food webs, on the other hand, empirical systems tend to have feasibility domains smaller than either their 
random counterparts or models conserving the directionality of interactions (cascade model). It is an open question 
which properties of real food webs are responsible for restricting the feasibility domain in this way. A possible 
candidate is the group structure observed in food webs |42j , corresponding to larger similarity of how certain species 
interact with the rest of the system than expected by chance, which in turn reduces the size of the feasibility domain 
(see Supplementary Information). 

These results parallel those for the distribution of the side lengths of the convex polyhedral cone delimiting the 
feasibility domain. The variance of side lengths for empirical structures is much higher than in random networks. 
This implies that even if the total size of the feasibility domain is large, it will have a distorted shape that is very 
stretched along some directions and shortened along others (Fig.[^. Consequently, it will be possible to find parameter 
perturbations of small magnitude that will drive the system outside its feasibility domain |43j . 

We have shown that each side of the feasibility domain corresponds to a pair of species, with the length determined 
by how similarly the two species interact with the rest of the system. As two species interact more and more similarly 
(i.e., have a larger niche overlap), the corresponding side becomes shorter and shorter, which in turns means greater 
sensitivity to parameter perturbations. Consistently with earlier results mi, this fact establishes a relationship 
between niche overlap and the range of conditions that lead to coexistence: greater niche overlap means a more 
restricted parameter range allowing for coexistence, irrespective of the details of the interactions. 

These differences between the size and shape of the feasibility domain shed light on the contrasting results obtained 
in the past on the effect of network structure on persistence [5DH5^ 15^ HOI HO] . Most of these studies rely on 
numerical integration, and therefore strongly depend on initial conditions. Given the difference in the shape of the 
feasibility domains of random and empirical networks, different initial conditions and their perturbations could result 
in markedly different outcomes: the feasibility domain could appear to be large or small depending on the direction 
in which perturbations are made. 


Methods 

Disentangling stability and feasibility 

From Eq. a feasible fixed point, if it exists, is given by the solution of 

s 

i=i 

where the asterisk denotes equilibrium values. A fixed point is locally asymptotically stable if all eigenvalues of the 
community matrix 

Mij = Uj Aij (4) 

have negative real parts. As discussed in the Supplementary Information, if A is diagonally stable or Volterra 
dissipative (i.e., there exists a positive diagonal matrix D such that DA + A^D is stable), then a feasible fixed point 


is globally stable in K+. 

A general characterization of diagonally stable matrices is unknown for S' > 3 [22]. There exist algorithms [44] 
that reduce the problem of determining if a S x S matrices is diagonally stable into two simultaneous problems of 
(S — 1) X (S — 1) matrices. While this method can be efficiently used to determine the diagonal stability of 4 x 4 
matrices, it becomes computationally intractable for large S. 

A matrix A is negative definite if 

XiA^jXj ^ 0 , (5) 

3 

for any non-zero vector x. A necessary and sufficient condition for a real matrix A to be negative definite is that 
all the eigenvalues of A + are negative [13]. A negative definite matrix is also diagonally stable, as the condition 
for diagonal stability holds with D being the identity matrix. Since it is extremely simple to verify this condition 
and it has been characterized for random matrices, we will study feasibility of negative definite matrices. In the 
Supplementary Information we show that with this choice we are excluding only a small region of the parameter 
space. 


Size of the feasibility domain 

The quantity S is the proportion of intrinsic growth rates leading to feasible equilibria. While a more rigourous 
definition is presented in the Supplementary Information, with a slight abuse of notation, S can be thought of as 

„ growth rate vectors leading to feasible equilibrium 

total # growth rate vectors 

The factor 2'® is an arbitrary choice that does not affect the results. It has been introduced to have 5 = 1 in absence 
of interspecific interactions (A^ = 0 if i 7 ^ j in Eq. and when all the species are self-regulated {An < 0 A i ^ j 
in Eq. [^. Given the geometrical properties of the feasibility domain, the proportion of feasible growth rates can be 
calculated considering only growth rate vectors of length one (Eig.j^and Supplementary Information), as this choice 
does not affect the value given by Eq. |S21| In the Supplementary Information we provide an integral formula for 
S (331 uni which makes both numerical and analytical calculations possible. 

Our method is still valid if some of the species are not self-regulated (i.e.. An = 0 for some i). In the Supplementary 
Information we explicitly discuss the properties of the feasibility domain of a community with consumer-resource 
interactions. In that case, S = 0 either when the diversity of consumers exceeds the diversity of resources or in the 
absence of interspecific interactions. Since consumers are regulated by their resources, they cannot survive in their 
absence and should therefore be characterized by negative intrinsic growth rates. We observe indeed that a necessary 
condition for an intrinsic growth rate vector to be contained in the feasibility domain, is to have negative values for 
the components corresponding to consumers. 
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Random matrices and moments 


El, E 2 , and Ec are moments of the random distribution for the off-diagonal elements of the interaction matrix, 
and are simply and directly related to the interaction strengths. They can be calculated as 

“ s(s-1) S ’ P) 

F - ^ A A 

^ ~ S(S - 1)E^ ^ El' 

For random networks with connectance C, these expressions reduce to [26] 

El =C^Ji, 

El = C{1 - C)fi^ + Ca^ , (8) 

F = + (1 - 

" + (1 _ 

where fi is the mean of the interaction strengths, cr is their variance, and p is the average pairwise correlation between 
the interaction coefficients of species pairs [26]. 


Universality of the size of the feasibility domain 

The size of the feasibility domain should, at least in principle, depend on all the entries of the interaction matrix. 
When these elements are drawn from a distribution, the size S of the feasibility domain is then expected to depend on 
all the moments of that distribution. As S increases, the dependence of S on some of those moments and parameters 
might become less and less important. S is universal if, in the limit of large S, it depends only on a few properties of 
the interaction matrix (i.e., on just the first few moments of the distribution). 

Specifically, for each unique pair of species (z, j), we set Aij = 0 with probability 1 — C and assign a random pair 
of interaction strengths {Mij, Mji) = {x,y) with probability C. The pair {x,y) is drawn from a bivariate distribution 
with given mean p,, variance a, and correlation p between x and y [26] . By considering different bivariate distributions, 
we can analyze the effect of different sign patterns (e.g., only (+, —) or (+, +) interactions) and different marginal 
distributions (e.g., drawing elements from a uniform or a lognormal distribution). 

Non-universality of 2 would mean that it depends on all the fine detaiis of the parameterization: 

2 = / (5', p, cr, p, C, sign pattern,...) , (9) 

where /(•) is an arbitrary function. The dependence on /i, u, and p can, without loss of generality, be expressed in 
terms of Ei, E 2 , and Ep- 

2 = p (S', Ai, A 2 , i?c, C”, sign pattern,...) . (10) 

However, if 2 is universal, then for large S, it is possible to express it as a function of Ei, £’ 2 , and E^ only: 


h{S,Ei,E2,E^) . 


( 11 ) 
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To verify this conjecture, we calculated S for matrices with the same values of Ei, £' 2 , and that differed for the 
values of the other parameters. As extensively shown in the Supplementary Information, S is uniquely determined 
by S, El, £ 2 , and £c (Eq. [^. 


Parameterization of mutualistic networks 


The 89 mutualistic networks (59 pollination networks and 30 seed-dispersal networks) were obtained from the Web 
of Life dataset (www.web-of-life.es), where references to the original works can be found. Empirical networks are 
encoded in terms of adjacency matrices L: Lij = 1 if species j interact with species i and 0 otherwise. When the 
original network was not fully connected, we considered the largest connected component. 

In the case of mutualistic networks, the adjacency matrix L is bipartite, i.e., it has the structure 


L = 




( 12 ) 


where Li, is a 5^ x Sp matrix {Sa and Sp being the number of animals and plants, respectively). The adjacency 
matrix contains information only about the interactions between animals and plants, but not about competition 
within plants or animals. 

We parameterized the interaction matrix in the following way: 


/ Lbo \ 

Ll o ) 


(13) 


where the symbol o indicates the Hadamard or entrywise product (i.e., {Ao B)ij = AijBij), while W^, 
and are all random matrices. and are square matrices of dimension Sa x Sa and Sp x Sp, while 
and are rectangular matrices of size Sax Sp and Sp x Sa- The diagonal elements and are set to — 1, 

while the pairs and {W^,W^) are drawn from a bivariate normal distribution with mean /i_, variance 

o'p = and correlation pa^. Since these two matrices represent competitive interactions, p- < 0. The pairs 
, Wjl^) were extracted from a bivariate normal distribution with mean pp, variance at = cp]^, and correlation 
pat, where pp > 0. For each network and parametrization we computed the size of the feasibility domain S. 

We considered different values of p-, pp, c, and p. Their values cannot be chosen arbitrarily, since A must be 
negative definite. Eor a choice of c, p, and a ratio p-/pp, the largest eigenvalue of [A + A^)/2 is linear in pp (as an 
arbitrary pp can be obtained by multiplying A by pp and then shifting the diagonal). Given the values of p-/pp, 
c, and p, one can therefore determine Pmax, the maximum value of pp still leading to a negative definite A (i.e., the 
value of Pp such that the largest eigenvalue of {A + A^)/2 is equal to 0). Fig. [^was obtained by considering more 
than 1000 parameterizations. Both the ratio p-/pp and the coefficient of variation c could assume the values 0.5 
or 2, while the correlation p assumed values from the set {—0.9, 0.5,0, 0.5,0.9}. The value of pp was set equal to 
0.25/irnax and 0.75/rniax- 
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Parameterization of food webs 

In the case of food webs the adjacency matrix L is not symmetric, Lij = 1 indicating that species j consumes 
species i. We removed all cannibalistic loops. Since Lij and Lji are never simultaneously equal to one (there are no 
loops of length two), we parameterized the offdiagonal entries of A as 

= W+L,, + , (14) 

while the diagonal was fixed at —1. Both and W~ are random matrices, where the pairs {W^, ) drawn 

from a bivariate normal distribution with marginal means (/r+,^_) and correlation matrix 

cul pcfil 

2 2 
pc^_ C}1_ 

We considered considering different values of /r_, c, and p. As explained above, given the values of p-lp,+ , c, 
and p, one can determine Pmax, the maximum value of still corresponding to a negative definite A. Fig. was 
obtained by considering more that 350 parameterizations. Both the ratio p_//i+ and the coefficient of variation c 
could assume the values 0.5 or 2, while the correlation p assumed either the value —0.5 or 0.5. The value of p,+ was 
set either to 0.25pmax or 0.75/i„iax- 
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FIG. 1: Geometrical properties of feasibility. The panels show the size and shape of the feasibility domain for three 
interaction matrices, each defining the interactions between three populations. If r corresponds to a feasible equilibrium, so 
does cr for any positive c; one can therefore study the feasibility domain on the surface of a sphere [25] (Supplementary 
Information). The gray sphere represents the S = 3-dimensional space of growth rates, while the colored part corresponds to 
the combination of growth rates leading to stable coexistence. The area (or volume for higher-dimensional systems) of the 
colored part is measured by H. Larger values of H correspond to a higher fraction of growth rate combinations leading to 
coexistence: the red interaction matrix is therefore more robust against perturbations of r than the green one. The size of 
this region (i.e., the value of H) does not capture all the properties relevant for coexistence. The red and blue systems have 
the same H, but the two regions-despite having the same area-have very different shapes, summarized in the bottom-right 
panel, where we show the length of each side for the red and blue systems. In the red system, the three sides have about the 
same length, and thus moving from the center in any direction will have about the same effect. In the blue system however, 
one side is much shorter than the other two, implying that even small perturbations falling along this direction may drive the 
system outside the feasibility domain. One of our main results is that, roughly speaking, if the red system corresponds to the 
random case, then the green one to food webs (having the same heterogeneity in side lengths as the random case but with a 
smaller H overall), and the blue one to empirical mutualistic networks (H rougly the same as in the random case but with the 

heterogeneity in side lengths much greater). 
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FIG. 2: Feasibility domain in random and empirical webs. The top two panels show H, the size of the domain of 
growth rates leading to coexistence, in the case of random networks. The left panel shows the dependence of H on Fi = C/r 
(where C is the connectance and /r is the mean interaction strength), and the number of species S. The right panel shows the 
match between our analytical prediction (Eq. [^and Supplementary Information) and the numerical value of H. The bottom 
panels show a comparison between H computed for empirical webs (89 mutualistic networks on the right, and 15 network was 
parameterized with different distributions of interaction strengths (Methods). Mutualistic networks have values of H 
comparable with random networks with similar interactions {B? = 0.98), indicating that their structure has little effect on the 
size of the feasibility domain. Food webs have lower values of H than their random counterparts (B? — 0.80). Empirical 
networks were parameterized extracting interaction strengths from a bivariate normal distribution with different means, 

variances, and correlations (Supplementary Information). 
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FIG. 3: Distribution of side lengths in random, structnred, and empirical networks. Left panels show the mean 
and the standard deviation of cos{ri), where rj is the side length. Analytical predictions for the first two moments of cos{ri) 
(Supplementary Information) perfectly match the numerical simulations. The two panels on the right show the standard 
deviation of cos{rj) for mutualistic and food webs compared to the expectations for the randomized cases. Both trophic and 
mutualistic interactions show larger fluctuations of side lengths, suggesting the existence of perturbation directions to which 
the system is more sensitive than to others. This effect is particularly pronounced and relevant for mutualistic networks. 
While mutualistic and random networks have a similar feasibility domain size H, this result implies that the response of 
mutualistic networks to perturbations is in fact more heterogeneous than those of their random counterparts. 
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SI. COMMUNITY DYNAMICS, FEASIBILITY, AND STABILITY 


We consider an ecological community composed of S populations, whose dynamics is described by the following 
equations: 


dn^ 

dt 


i=i 


V'V j ’ 


(SI) 


where is the population abundance of species i and is its intrinsic growth rate, and Aij is the effect of a unit change 
in species j’s density on species i’s per capita growth rate. For notational convenience, we collect the coefficients Aij 
into the interaction matrix A, and nt and into the vectors n and r, respectively. 

In principle, the interaction matrix A may depend on n. We discuss this more general case in section |S14[ In 


the following, we consider the simpler case of A being independent of n; then, equation (SI) is a general system of 
Lotka-Volterra population equations. 

A vector n* is a fixed point (equilibrium) if 

s 

0 = + (* = I,2,...,5) . (S2) 

A fixed point is feasible if n* > 0 for alH. A feasible fixed point (if it exists) is then a solution to the equation 

s 


i=i 


(S3) 


and therefore, assuming A is invertible, 


= -E(^ ■ 

i=i 


(S4) 


A fixed point n* is locally stable if the system returns to it following any sufficiently small perturbation of the 


population abundances. Introducing rii = n* + 6ni in equation SI and assuming that Srii is small, we obtain, by 
expanding around Sn^ = 0, 

s 


dSrii 

dt 




(S5) 


i=i 


where is the {i,j)th entry of the Jacobian evaluated at the fixed point (also called the community matrix), which, 
in the case of equation [ST| reduces to 

s 

■ 


Mij = n*Aij = - [ ^(A '^)ikrk ) A, 


(S6) 


\k^l 


Substituting into equation S5, we get 


dSrii 

dt 


= - X! (E^^ A,jSnj 

j=i \k=i J 


(S7) 


There are two possible scenarios for the dynamics of equation |S5| If all eigenvalues of M have negative real parts, 
then the perturbation 6n decays exponentially to zero and n* is locally stable. If at least one eigenvalue of M has a 
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positive real part, then there exists an infinitesimal perturbation such that the system does not return to equilibrium. 
If we order the eigenvalues Xi of M according to their real parts, i.e., 3fJ(Ai) > 3fJ(A2) > • • • > 5i(As), then stability 
depends exclusively on 3?(Ai): if it is negative, n* is dynamically locally stable; otherwise, it is unstable pS] , 

A fixed point is globally stable if it is the final outcome of the dynamics from any initial condition involving strictly 
positive population abundances. 


S2. DISENTANGLING STABILITY AND FEASIBILITY 


As we can see from equations S4 and |S7[ both feasibility and stability depend on both r and A and, at least in 
principle, a fixed point can be stable or unstable, independently of the fact that it is feasible or not. 

We want to study the proportion of conditions (i.e., the number of combinations of the growth rates r out of all 
possible combinations) leading to coexistence, i.e., leading to stable and feasible equilibria. Therefore in principle we 
should, for a fixed matrix A, look for growth rates r that satisfy both stability and feasibility. In probabilistic terms, 
we want to measure the likelihood that a random combination of the intrinsic growth rates corresponds to a stable 
and feasible solution. 

In the case of equation |S1[ it is possible to disentangle feasibility and stability by applying a mild condition on the 
interaction matrix A. To this end, we introduce some terminology |47[ section 2.1.2]: 

• Stability. A real matrix S is stable if all its eigenvalues have negative real parts. 

• iD-stability. A real matrix S is D-stable ii D B is stable for any diagonal matrix D with strictly positive 
diagonal entries. 

• Diagonal stability. A real matrix B is diagonally stable if there exists a positive diagonal matrix D such that 
D B + B^ D is stable (where B^ is the transpose oi B). 


We also consider 


• Negative definiteness (in a generalized sense). A real matrix B is negative definite if J2ij ^iBijXj < 0 for 
any non-zero vector x |23j . 

These properties are closely related to each other mm- 

Negative definiteness Diagonal stability D-stability Stability (S8) 


• Negative definiteness => Diagonal stability. A matrix B is negative definite if and only if all the 
eigenvalues of B -f B^ are negative [23] . If this condition hold, then the positive diagonal matrix satisfying the 
definition of diagonal stability is simply the identity matrix. 

• Diagonal stability => D-stability. See the book by Kaszkurewicz & Bhaya for the proof jU] lemma 2.1.4]. 

• D-stability Stability. This follows from the definition of D-stability when D is the identity matrix. 

In the case of equation |S1[ those conditions applied to the matrix A are related to the stability of the system. One 
can use the definition of the community matrix (equation |S6|) to show that D-stability of A implies the local 
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asymptotic stability of any feasible fixed point. This is because the community matrix with entries = n*Aij 
can be written as N A, where N is the diagonal matrix with Na = n*. If the fixed point is feasible and A is D-stable, 
then local asymptotic stability is guaranteed. Moreover it is possible to show [SI in] that diagonal stability of A 
global stability. 

Thus, we have a condition on A that makes it possible to disentangle the problems of stability and feasibility: 
A is negative definite => global stability of the feasible fixed point [50]. Therefore, if we assume A is 
negative definite, then feasibility of the equilibrium is sufficient to guarantee its global stability as well, i.e., feasibility 
guarantees globally stable coexistence. Consistently with this, it is known that the largest eigenvalue of {A + A^)I2 
is always larger than or equal to the real part of A’s leading eigenvalue [53], i.e. negative definiteness implies stability. 
While this was indeed observed before, it is important to underline that, in the case of ref. (53], this property was 
considered on the community matrix M (which also depends on the fixed point’s position in phase space) and not on 
the interaction matrix A. 

Since we are interested in studying how interactions (i.e., the matrix A) determine coexistence, and which properties 
of the former determine the latter, we will restrict our analysis to negative definite matrices A and focus only on the 
problem of feasibility. This condition has the advantage of being analytically computable for large random matrices 


(see section S5 A). 


S3. GEOMETRICAL PROPERTIES OF THE FEASIBILITY DOMAIN 


In section [S2] we showed how to separate feasibility and stability, i.e., we have a sufficient condition on the interaction 
matrix that guarantees (global) stability of the feasible fixed point. The problem of determining the size of the 
coexistence domain is therefore reduced to that of determining the size of the feasibility domain. The ecological 
interpretation of this volume is the proportion of different conditions leading to feasible equilibria out of all possible 
conditions. The larger this volume is, the higher the probability that the system is able to sustain biodiversity. In 
terms of equation |SH we want to quantify the proportion of growth rate vectors r corresponding to a feasible hxed 
point. 

This geometrical approach was pioneered in |25] where the space of feasible solution was studied for dissipative 


systems, and the size of that domain was computed in the case S = 3 (see section S13|. 

At this point, it is important to observe that if a vector r corresponds to a feasible solution, then cr, c being an 
arbitrary positive constant, also corresponds to a feasible solution. This is because the equilibrium solution n* is 


given by equation S4 which is linear in r^. Therefore, the equilibrium corresponding to cr^ is simply cn*, and since c 


is positive, cn* is also feasible. 

This fact implies that, given a large number of growth rate vectors r, the expected proportion of vectors correspond¬ 
ing to a feasible fixed point is independent of r’s norm. In other words, r is feasible if and only if r/||r|| is feasible, 
where ||r|| = is the Euclidean norm of r. The proportion of feasible growth rates among all possible ones is 

therefore equal to the proportion of feasible growth rates calculated using only growth rate vectors with ||r|| = 1; i.e., 
those lying on the unit sphere. 

Before proceeding with the mathematical definition of the size of the feasibility domain, we discuss the geometrical 
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interpretation of equation S4 


From this equation, the feasibility condition reads 

s 


< 0 . (S9) 

This equation defines a convex polyhedral cone in the S'-dimensional space of growth rates. A convex polyhedral 
cone [IH] is a subset of M'® whose elements x can be written as positive linear combinations of Nq different S- 
dimensional vectors called the generators of the cone: 


Ng 

X = ^g'^Xk , 


(SIO) 


where the Afc are arbitrary positive constants. Due to this arbitrariness, if is a generator of a given convex 
polyhedral cone, then also cg^ (where we rescale just the /cth generator with the positive constant c, leaving the 
others unchanged) will be a generator of the same cone In the case of equation S3 each and every growth rate 
vector belonging to the feasibility domain can be written as 

s 

r^ = -'^Aiknl, (Sll) 


where, by definition, is feasible and therefore a positive constant. One can easily see that this equation corresponds 


to equation SIO where the number of generators Nq is equal to S and the ith component of the vector g^ is proportional 
to —Aik- As the lengths of the generators can be set to any positive value, we will normalize them to one, i.e.. 


gf(A) = . (S12) 

The generators completely define the feasibility domain in the space of growth rates. A growth rate vector corresponds 
to a feasible equilibrium if and only if it can be written as a linear combination of the generators with positive 
coefficients. Biologically the generators correspond to the growth rate vectors that bound the coexistence domain. 
They correspond to nonfeasible equilibria with just one species with positive abundance (and all the others with zero 
abundance), such that there exist arbitrarily small perturbations of the growth rate vector that make the equilibrium 
feasible. 

The set of all the growth rate vectors leading to a feasible equilibrium is therefore a convex polyhedral cone, defined 

by 


s 

K{A) = {rGR^\J2{A-%r,<0} . 


(S13) 


Equivalently, it can be defined in terms of the generators: 


s 

K{A) = {rGR^\3Xi,X2,...,Xk>0,r = Y,gHA)Xk} , 


(S14) 


k=l 

where the generators g^{A) are defined in equation 
out in the case of S' = 3. 

This geometrical definition and characterization of the feasibility domain allows us to identify classes of matrices 
having the exact same feasibility domain: they are simply matrices having the same set of generators. In particular. 


SI2 


In section 


SI3 we show explicitly how these concepts pan 
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there are two basic transformations of the matrix A (and their combinations) that leave the set of generators un¬ 
changed: permutations and positive rescaling. A square matrix P is a permutation matrix if each row and column 
has one and only one nonzero entry and the value of that entry is equal to one. A positive rescaling is performed by 
a positive diagonal matrix D. The set of generators of A is the same as those of P A and D A. This can be seen by 
observing that a permutation of the rows just changes the order of the generators but not the generators themselves. 
In the same way, a generator with the same direction but different length generates the same cone, and so any positive 
constant that rescales a row of the matrix leaves the feasibility domain unchanged. It is important to note however 
that these two transformations do not leave the properties of the matrix A unchanged: both exchanging rows of a 
matrix and rescaling rows by different constants will in general change the structure of the matrix. 

Using this geometrical framework, one can easily identify the center of the feasibility domain (also known as 
structural vector i)- There are several possible ways to define the center of a hypervolume and, without additional 
assumptions, all the definitions are different. One natural choice is the barycenter (“center of mass”) of the domain of 
feasible intrinsic growth rates. Any plane passing through the barycenter divides the volume into two subvolumes of 
equal size. The barycenter is equivalent to the center of mass of the volume (in the case of constant density). Then, 
the vector pointing from the origin to the barycenter is given by 


x'’ = 


y 


(S15) 


JK{A)nSs 

where 0 is the intersection of two sets, and Sg = {r e IR'^|||r|| = 1} represents the surface of the S'-dimensional unit 
sphere. The variable y is therefore integrated over the feasibility domain restricted to the unit sphere’s surface. All 
points in the feasibility domain are positive linear combinations of the generators, i.e.. 




(S16) 


where the are positive constants. The fact that we consider only the points lying on the unit sphere, i.e., ||y|| = I, 

as 


can be expressed as a constraint on A (the vector of As). Thus, we can write equation SI5 

= [ d^Ag(A)^AV , 


X 


(S17) 


where q is an appropriate distribution, introduced to take into account three different constraints: all the components 
of A must be positive; the vector must lie on the unit sphere; and those vectors must be sampled uniformly 

on the feasibility domain. One can show that the distribution q{X) has the following form 


q{X) oc exp - ^ A*(g* • g^)X^ ©(A'”) 


* J / k 

where the proportionality constant is given by the normalization. Therefore, by defining. 


we obtain 


J d^Xq{X) A'==: (A'^) , 




(S18) 


(S19) 


(S20) 
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S4. DEFINITION AND CALCULATION OF H 


As explained in section the proportion of feasible growth rates can be calculated considering only growth rate 
vectors of length one, i.e., ||r|| = 1. This proportion can be interpreted as the volume of the intersection of a convex 
cone and the surface of a sphere. Equivalently, it is the solid angle of the convex polyhedral cone [iHl fiH] , 

We define the quantity 5 as 

„ # growth rate vectors corresponding to a feasible fixed point (S21) 

total # growth rate vectors 

The factor 2^ that appears in this equation is an arbitrary choice, and it has been introduced to have S = 1 when 


species are not interacting = 0 if i ^ j). In this case equation SI reduces to S independent logistic equations with 


equilibrium densities n* = —rifAu. Taking each An to be negative (otherwise each species would have an unstoppable 
positive feedback on itself), this equilibrium is feasible if and only if each is positive. For a single species then, 
the probability of randomly drawing a feasible (i.e., positive) growth rate out of all possible growth rates is one half. 
For two species, both growth rates must have the correct sign to have the two species with positive abundance, and 
therefore the proportion of growth rate vectors satisfying this condition is 1/4. For S species the combinations of the 
growth rates leading to a feasible fixed point is 2~^. S, defined as in equation S21 is therefore equal to one when 
species do not interact. 


In terms of geometrical properties and the convex polyhedral cone, S can be defined as 

c.vois_i(Ar(A) n Ss) 


= 2 " 


V0l5_i(Ss) 


(S22) 


where K{A) is defined in equation S13 Sg is the unit sphere in K*, while vols(-) means volume in S dimensions. 


This definition is equivalent to the one in equation S21 [45l 146] . 


These two equivalent definitions can be expressed in terms of an integral in the space of the growth rate vectors: 

oS - S 


V0ls_i(Ss) 7 rs 


d‘®r 2||r|| 


|rf - 1) ]j0(n*(r)) 


(S23) 


where vols_i(Ss) is the volume of the unit sphere’s surface in S dimensions, ©(•) is the Heaviside function (equal to 1 
is the argument is positive and to zero otherwise), and S{-) is the Dirac delta function. In this expression, we integrate 
over the surface of the S'-dimensional unit sphere. The integral of a function f{x) on the unit sphere is given by 


[ d'^a: f{x) = / d‘ 

JSs 4ms 


X 2\\x 


- 1 ) fix) 


(S24) 


where the term 5(||a;||^ — 1) that appears in the integration constrains x on the surface of the unit sphere, and the 
factor 2||a;|| is the derivative of the delta function’s argument, which is needed because the Dirac delta is nonlinear in 
||r||. The factor vols_i(Ss), the surface of sphere in S dimensions, can be obtained by setting f{x) = 1: 

27r«/2 


vols_i(Ss) = y 


X 2\\x 


- 1 ) = 


r(5/2) 


(S25) 


where r(-) is the Gamma function. Finally, the term ('^)) ™ equation S23 expresses the constraint of all 


n* having to be positive: this product is equal to 1 if the equilibrium n*{r) is feasible and zero otherwise. The 


equilibrium n*(r) is a function of r via equation S4 














23 


Equation S23 defines S as the volume of the domain of growth rates leading to feasible solutions. Using the results 
of section [S^ we know that if the interaction matrix A is negative definite then a feasible fixed point is globally stable. 
In this case S is the volume of the domain of intrinsic growth rates leading to feasible and (globally) stable solutions. 

Unfortunately, direct numerical computation of 2 is inefficient when the number of species S is large. To evaluate 
the integral in equation |S23[ e.g., via Monte Carlo integration, we should draw intrinsic growth rates at random and 
count how many of them, out of the total, lead to a feasible equilibrium. In order to have a reliable estimate of this 
proportion, we should sample the space in such a way that the number of feasible growth rates found is large. This 
goal requires an exponentially increasing sampling effort as S increases. In this section we provide an alternative, 
much faster and reliable, way of estimating 2. 

The equilibrium solution and the growth rates are linearly related via ri = — (equation 


S3). Our 


strategy is to use this to perform a change of variables in equation S23 and integrate over n* instead of r. Since 
A is negative definite (and thus stable and not singular), it is invertible, and so it is always possible to perform this 
change of variables. Note that, more generally, the change of variables can be performed if A is nonsingular (i.e., 
det(yl) = 0). We then obtain 


2^ r(5’/2) |det(A)| 


f d^n* 26ly^n*Ak,A 


kjn* - 1 


nQ«) ’ 


(S26) 


where | det(A)| is the determinant of A, which is also the Jacobian of the change of variables. After the change of 
variables, the integration is now performed over the feasible equilibrium points and so the condition of feasibility is 
automatically implemented. 

It is still difficult to evaluate the previous expression numerically, because of the constraint that appears in the 
delta function. We can further simplify it by introducing polar coordinates. In particular, we write the vector n as 
n = nu, where n = HnH and m is a vector of unit length. We can perform a new change of variables, passing from n 
to n and u. Specifically, for any function f{n), we can write 

p poo p poo p 

/ d^n f{n) = dn / d^u 2J(||m|P — l)f{nu) = dn / d^u f{nu) . (S27) 

Jrs Jo Jrs Jo Jss 


Using this expression in equation |S26[ we obtain 
„ 2® r(5/2) det(A) 


2'k^/'^ 


dn 


s-i 


f d^u 26 i V? ^ UiGijUj ~ 11 

\ i,j ) i=l 




(S28) 


where we used the fact that ©(n^) = 0(ui) (since rii = nui, and n is positive by definition), and we have introduced 
the matrix Gij = A^iAkj. We can now perform the integration over n, obtaining 


[ dn n®-i 2 (5 I 

UiGijUj - 1 

Jo ' 



-S/2 


(S29) 


dn n'^ ^ 2 din — 


j UiGijUj 


2 n ^i j UiGijUj 






and therefore the integral of equation |S23| finally reads 

„ _ 2^ r(5/2) ^det(G) 


-s /2 


2'it^G 


d^u ]^ 0 ('u^) y^^UjG.jUj 


(S30) 


2 = 1 
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where we have used the fact that det(G) = det(A^A) = det(A)^. In terms of the interaction matrix, the equation 
reads 


_ _ 2^ r(5'/2) |det(A)| 


-S/2 


d'^M Q{ui) I ^ UiAkiAkjUj 


(S31) 


Equation S30 shows explicitly the role of the generators. The matrix G can indeed be rewritten as 


G 


ik 


= X! 9j9jCiCk = c^Ckg'- ■ g’^ 


(S32) 


where are the generators of the convex cone defined in equation 


S12 


and Ci are arbitrary positive constants. Their 
presence, which can be seen as a change of the normalization of the vectors does not affect the form of equation S30 


and its dependence on G (see section S31. This property can be checked explicitly from equation S30 by introducing 
an explicit dependence on Ci and showing that S is independent of their values. 

Unfortunately, the integral in equation |S30| cannot be computed analytically. As mentioned before, when the 


integral is written in the form of equation S23 it is impractical to evaluate it numerically, since it would require an 
exponentially increasing sampling to get a reasonable precision. Fortunately, this is not the case when the integral is 
written as in equations |S30| and |S31[ The main difference is that, after changing variables, we are directly sampling 
the space of feasible solutions, without losing computational time in randomly exploring the space of intrinsic growth 
rates looking for feasible solutions. 

To evaluate the integral, we use the usual approach of Monte Carlo algorithms. In particular, it is possible to write 
the integral as an average over random points: 

T am ^ S 


1 

f 




-s /2 


ns/2) 

2Tr^G 


d^u Pj0(ui) 2(5 (||m|P - I) (^u,GijU^ 


-s /2 


_ (S33) 

a —1 i,j 2=1 ij 

when r —>■ oo. In this expression are independently drawn random vectors uniformly distributed on the unit sphere 
and with only positive components. These two conditions are introduced to satisfy the constraints Ilf=i 0(^i) 

2(5(11 u|p — 1) that appear in the integral. T is the sample size, and the average on the left hand side of equation S33 
converges to the right hand side in the large T limit. 

One always has a finite sample size T, used to approximate the integral. It is therefore important to have an 


estimate of the error made due to T < oo. Since the left hand side of equation S33 is an average of a function over 
random vectors, this error can be estimated by simply using the variance of the function’s values. In particular, the 
error ctmc is defined as 




I 

f 


E 





(S34) 


The numerical simulation presented in the work where obtained were obtained with different sampling effort T. 
Instead of fixing T a priori, we determined a precision goal, that we measured in terms of the relative error ctmc/S. 
We ran the simulations until (Tmc/ 5 < 0.05. In order to avoid artificially small samples and to have enough statistical 
power not to undershoot to much (TmCi ran 10 x S' Monte Carlo steps before checking the condition for the first 
time. 
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S5. STABILITY, NEGATIVE DEFINITENESS, AND FEASIBILITY IN RANDOM MATRICES 

Random matrices are a useful tool in ecology, and have been studied since May’s seminal paper m- Mostly, they 
have been used to model the community matrix [nmi]. In the context of this work, we use random matrices to 
model interaction matrices A. We consider random matrices constructed in the following way: 

• An = —d where d is a positive constant. 

• Each pair is set equal to a pair of random variables drawn from a joint distribution with probability 

density function q{x,y). 

• The random variables are exchangeable—i.e., the probability distribution function is symmetric in its arguments: 
q{x,y) = q{y,x) —and all the moments are finite. 

We show that the three most important quantities for our problem are the moments 

Ei= dx dy xq{x,y) = dx dy yq{x,y) , (S35) 


Eo = 


^ j dx dy {x- EiYq{x, y) = ^J dx dy {y - Ei)'^q{x, y) 


(S36) 


Ec= ^ J dx dy {x- Ei){y - Ei)q{x,y) . 


In the limit of large S, they can be computed as proper sample means of A’s entries: 




(S37) 


(S38) 


E2 


\ S{S-1) 




i=l j^i 


(S39) 


Er = 


El \SiS-l) 




i=l j^i 


The parameterization used by May m would correspond to 


(S40) 


qMa,yix, y) = ((1 - C)6{x) + Cp[x)^ ((I - C)5{y) + Cp(y)^ , (S4I) 

where d(-) is the Dirac delta function and p{x) is an arbitrary distribution with mean zero and variance a^. The 
connectance C sets the probability that each entry is equal to zero (with probability I — C) or randomly drawn from 
the probability distribution p{x) with probability C. In this case Ei = Ec = 0, while El = Ca^. 

In the following, we summarize known results on the spectra, negative definiteness conditions, and properties of 5 


for these matrices. 
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A. Known results on the spectra of random matrices 

Under the assumptions of the previous section, the eigenvalues of A in the limit of large S are uniformly distributed 
in an ellipse in the complex plane. If 7 ^ 0 there is always an eigenvalue Am whose value is approximately 

Am « -d + , (S42) 

independently of the rest of the eigenvalue distribution. The ellipse is centered at —d — Ei, its axes are aligned with 
the real and imaginary axes, and their lengths are 

a = '/SE2{l + E^) (S43) 

and 

h = y/SE2{l-E^) . (S44) 

If Am = 0, the eigenvalue with the largest real part(s) is approximated by the rightmost point of the ellipse. The 
system is stable if its real part is negative. In the most general case, this condition is equivalent to 

-d + max|s'Ui,-Ui+ v^£;2(1 + Uc)} <0 . (S45) 

In section |S2| we introduced the concept of negative definiteness. In particular, we showed that when the matrix 
is negative definite then it is possible to disentangle stability and feasibility. The matrix is negative definite if the 
eigenvalues of A + are all negative. This condition reads [23] 

-d + max|5'Ui,-Ui + y/25'(l + T;e)U 2 | <0 . (S46) 

Figurej^shows the values of parameters leading to the possible combinations of stability and negative definiteness in 
random matrices for the case Ei = 0. Since we imposed that A is negative definite, the region of parameters we explore 
is the one above the negative definiteness line. One can see that in this way we are missing some parameterizations, 
corresponding to those that lead to a stable but not negative definite matrices. From equations |S45| and |S46| one can 
see that the case Ui < 0 is very similar to the case Ei = 0. More interestingly, for Ei > 0, the conditions for stability 
and negative definiteness converge in the large S limit, implying that we are considering all the possible cases. 

What is remarkable in these conditions and in the distribution of eigenvalues is that they are universal [2H [5DH5^ . 
Universality means that they depend only on S', Ui, U 2 , and Ec (and d, but via a trivial dependence). The spectrum 
of eigenvalues does not depend on the detailed form of the distribution q(x,y). 

For instance, consider the case q{x,y) = p{x)p{y), where the upper and lower triangular entries A^- and Aji are 
independent random variables. In this case Ec = 0 and Ei and E 2 are the mean and standard deviation of the 
distribution p{x). The distribution of eigenvalues and the conditions for stability and negative definiteness are the 
same for any probability distribution p{x) as long as their mean Ei and standard deviation E 2 are the same (provided 
some mild conditions on higher moments hold). For instance, a Lognormal distribution, a Gaussian distribution and 
an exponential distribution, having same mean and standard deviation, produce the same eigenvalue distribution, 
and therefore the same conditions for stability [53| . 

From an ecological perspective, one can consider different interaction matrices corresponding to different interaction 
types. The interaction type is given by the signs of the pairs {Aij, Aji): competitive interactions will have both entries 
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Supplementary Figure SI: Negative definiteness and stability for random matrices in the case Ei = 0. The red curve 


describes the condition for stability (equation S45l, while the blue curve corresponds to the negative definiteness condition 


(equation S46l. The region above the blue curve corresponds to matrices that are both stable and negative definite, while the 
region below the red curve corresponds to unstable and non-negative definite matrices. The parameterizations that may still 
lead to stable and feasible points but we are not considering are in the region between the two curves. The shape of this 
region does not change substantially if S and E 2 are changed or if Ei < 0. For Fi > 0 the not negative definite but stable 
region is always smaller and eventually disappears (i.e., the blue and the red curve become the same) when S is large enough. 


with a negative sign, while in trophic interactions the entries will have opposite sign. The interaction pairs 
for competitive interactions can for instance be obtained from the following distribution: 

qcomp{x,y) = {1 - C)6(x)6{y) + Ch_{x)h_{y) , (S47) 

where h_ is a probability distribution function with support on the negative axis (i.e., the random variables are 
always negative), and C is the connectance (a pair is different from zero with probability C). In the case of trophic 
interactions we could consider 

qtroph{x,y) = (1 - C)6{x)6{y) + ^p-(x)p+{y) + ^p+{x)p-{y) , (S48) 

where and p- are two probability distribution functions with positive and negative support, respectively. Suppose 
that the moments of h_, and are chosen in such a way that gcomp(a:, y) and (i'troph(a^, y) have the same values 
of El, E 2 , and Ec- The interaction matrices will still look very different in the two cases: one describes a foodweb 
and the other a competitive system. Despite this difference, the two will have the same stability properties. In other 
words, different interaction types influence the stability properties of the system only via Ei, E 2 and E^. 
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B. Universality of H 

In this section we show that, apart from their spectral distribution, S is also a universal quantity in large random 
matrices. That is, in the large S limit, its value does not depend on the entire distribution of the coefficients, but 
only on the three moments Ei, E 2 , and Ec- It is important to remark that this result applies to the large S limit: 
the sub-leading corrections depend in principle on all the moments. 

In order to show that 2 is universal, we parameterized random networks with different distributions and checked 
whether 2 depends only on Ei, E 2 , E^ and S, but not on other properties. To do this, we constructed several S x S 
matrices. Each individual matrix had its entries drawn from some fixed distribution, but the shape of the distribution 
was different across matrices. However, regardless of the distribution’s shape, their moments were fixed at Ei, E 2 , 
and Ec- We then checked whether these matrices led to the same value of 2. 

In our simulations we considered a distribution of the pairs [Aij,Aji) of the form 

y) = {l- C)5{x)5{y) + Cp{x, y) , (S49) 

where the connectance C is the probability that two species i and j interact. The probability distribution p{x,y) in 
equation |S49| depends on three parameters p, a, and p, which define the mean, variance, and correlation of the pairs 
drawn from p{x,y). Given the values of i?i, 1^2, and Ec, we can arbitrary choose C and tune p, a, and p to obtain 
any desired Ei, E 2 , and Ec- If S is universal, then different matrices built with different values of C, p, a, and p but 
the same values of Ei, E 2 , and Ec will lead to the same 2. 

We considered five parameterizations of the distribution p{x,y): 

• Random signs, normal distribution: 

p{x, y) = BN{x, y\p, a, p) . (S50) 

The distribution BN{x,y\p,a, p) is a bivariate normal distribution with marginal means equal to p, marginal 
variances equal to a^, and correlation equal to pa^. The pairs can in principle assume all possible combinations 
of signs. 

• Random signs, four corners: 

p{x, y) = 7 :S{x - p- a)S{y - p-a) + ^6{x - p + a)6{y - p + a) 

1-q I-g 

H- - p- <^)Hy - p + cr)-\ - - P + cr)S{y - p - a) . 

The pairs {x,y) can take on only four different, discrete values, potentially corresponding to all combinations 
on signs. The probability distribution depends on three parameters p and are means an variances of the 
distribution, while the correlation pa'^ can be obtained from p = 2q — 1. 

• (-I-,-k), Lognormal: 

p(x, y) = LBN{x, y\p, a, p) . (S52) 

The distribution LBN{x,y\p,cr, p) is a bivariate lognormal distribution with marginal means equal to p > 0, 
marginal variances equal to cr^, and correlation equal to The pairs can in principle assume only positive 
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signs. Note that not all values of p between —1 and 1 can be obtained when a Lognormal distribution is 
considered. 

• Lognormal: 


p{x, y) = LBN{-x, -y\ - p, a, p) . 


(S53) 


This distribution takes the values drawn from a bivariate lognormal distribution, times —1. It has marginal 
means equal to /r < 0, marginal variances equal to cr^, and correlation equal to pcr^. The pairs assume only 
negative signs. Note that not all values of p between —1 and 1 can be obtained when a Lognormal distribution 
is considered. 


• (+,—), Lognormal: 


Pi.x, y) =]^LN{x\pi, (1 + p)a)LN{-y\ - p 2 , (1 + p)cr) 

+ ^LN{y\pi, (1 + p)a)LN{-x\ - p 2 , (1 + p)cr) . 

The distribution LN{x\p, a) is Lognormal distribution with mean pi + p 2 (where pi > 0 and p 2 < 0), 
cr^, and correlation pa'^. The pairs assume only values with opposite signs (+, —) or (—, +). 


(S54) 

variance 


In ecological terms, the first two distributions correspond to a random community (where the signs of the interaction 
strength are random), the (+,+) case corresponds to a mutualistic community, (—,—) to a competitive community, 
while (+,—) corresponds to a food web. The mutualistic/competitive matrices can lead only to positive/negative 
means Ei , respectively, while the other settings can produce arbitrarily values of Ei. 

Figure shows the value of S and of the largest eigenvalue A for interaction matrices constructed with different 
connectances C and distributions, but with the same values of ifi, E 2 , and E(.. As seen from the figure, the values 
of S and A in any particular case match up precisely with the average values over several different realizations, 
demonstrating that these two quantities are indeed universal. 


S6. MEAN-FIELD APPROXIMATION OF H 


The goal of this section is to compute an approximation for S in the limit of large S. The volume is defined (see 
section S4| as 

-S/2 


_ _ 2^ r(S'/2) v/det(G) 


d'^M ]^0(uj) y^UjG. 


n ^3 


(S55) 


where the matrix G can be obtained from the generators of the polytope (see equations S12 and S32|, and therefore 
from the interaction matrix A. 


We can introduce a Gaussian function in equation |S55| using the fact that, for any positive constant c, 

2 /’°° 

nW)L 


r.-S/2 _ 


(S56) 
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Supplementary Figure S2: Universality of A and H in random matrices. The two left panels refer to the eigenvalue with 
the largest real part A of the interaction matrix A, while the right ones to the size H of the feasibility domain. We consider 
different values of the connectance (colors) and different distributions (shape), such that there were multiple combination of 
connectances and distributions having the same values of Ei, E2, and E^. We computed the averages (A) and (log(H)) over 
all realizations of the matrices having the same values of Ei, E2, and E^. If the value of A and H are universal, then they 
depend only on Ei, E2, and Ec, and therefore their values are equal to the mean: universality holds if A = (A) and 
log(H) = ( log(H)). The top panels show that these two quantities are equal and the bottom panels quantify their deviations. 

We know that A is universal, and since H has a similar behavior, we conclude that H is also universal. 


Introducing this Gaussian integral in equation S55 by letting c = j UiGijUj, we obtain 

S = -ydet(G) J dr J d^u exp |^ UiG. 


U.J 


which can be rewritten as 


0 s r ^ 

\/det(G)(^^^ J ^d^z exp(^-Y^ZiGijZj'^ , 

i — 1 i^j 


(S57) 


(S58) 
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where Zi = rui. We can rewrite this equation as 

s 




Q{zi) e 



(S59) 


where we used the fact that the diagonal entries of G, when expressed in terms of the normalized generators, are 
equal to one. 

The reader familiar with statistical mechanics will notice that equation |S59[ which can be written as 

d^z q{z) n 1 1 “ X! 


/KS 


(S60) 




has the form of a partition function. For instance one can recover the Ising model [54j with the choice q{z) = 
Hi S{zf = 1) or the spherical model [53] when q{z) = S{S — ^ ■ zf). The term ZiGijZj in particular plays the role of 
the interactions of the system. 

Integrals of the form |S60| are the most studied objects of statistical mechanics, and yet in most cases are not 
analytically solvable. There are, on the other hand, many techniques that can be used to obtain good approximations 
to S60 The most celebrated one is probably the mean-field approximation |54j and it is the one we are using in this 


section. In particular, the idea of the mean-field approximation is to replace the interactions of an entity (spins in 
the case of the Ising model or species in our case) with an average “effective” interaction. This reduces a many-body 
problem, where all interactions of spins or populations are coupled, into an effective one-body problem. 

If the system is large enough (in our case if S' —>■ oo), the mean-field approximation is know to be exact in the case of 
“fully connected” interactions. In terms of equation |S60[ this corresponds to a matrix G with the same constant in all 
its offdiagonal entries. The matrix G is constant when A has constant offdiagonal entries. We will consider therefore 


the case of A’s diagonal entries being equal to — I and its offdiagonal entries to a constant Ei. Using equation |S12 
the ith component of the fcth generator is then 

El 


for i ^ k, and 


9i = 


gl = 


I + (S - l)El 


I 


I + (S - l)El ■ 


(S6I) 


(S62) 


Using equation S32 we therefore obtain that the diagonal entries of G are equal to I, while the offdiagonal ones are 
constant and equal to 

-2Ei + {S- 2)El 


We define the constant /3 as 




/3 = 5 . 


I + (5 - l)El 

-2Ei + {S- 2)El 
l + {S-l)El 


and therefore we have Gu = I and Gij = (3/S for i ^ j. The determinant of G in this case turns out to be 


det(G) — ( I H- ^—P 




s-i 


(1 + / 3 )£ 


(S63) 


(S64) 


(S65) 
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where the last form holds for large S. In this case of constant interactions, we obtain, from equation S59 

s 


^ =^/det{G) ( ^ j I 0(zi) e exp ( -2*^ XI 


i=l 


3¥=i 


— \J det(G) ( —!= I I 

, V TT / Jrs 


d^z exp ( 


up to subleading terms in S. 


Equation |S66| can be written as 


where 


= \/det(G) ^/i(exp ^-^(X^i)^+ ^X- 


Zh ■■= d®2: [ e(zi) ) exp [ - X ^ X = 


Jrs 


/ \ i i 

s s 

dz eik{h/2) 


where erfc(-) is the complementary error function, defined as 

erfc(a;) = / dt 

V'^ Jx 

The average {■)h is defined as 


exp (-X^*^ - *x^0 ■ 




Using Jensen’s inequality in equation |S70| we have that 


^ =v/det(G) ( X ) ^/»(exp (-|(X^i)^+ ^X^* 


S'' 


> 


> A/det(G) ( X ) ZheiipU-^{^z,f + h^Zi 


(S66) 


(S67) 


(S68) 


(S69) 


(S70) 


(S71) 


In the following we will approximate the first expression with the second one. It is possible to prove that, in the large 
S limit, the second expression converges to the first one. 

Applying the mean-field approximation we neglect fluctuations of the variables, i.e. we have 


+ + ~ ^ {-Prn^ + hm) , 


(S72) 


where 


1 d 

m := {zi)h = -^^\og{Zh) . 


By introducing equation S72 in equation S7I we have 


s « \/det(G)Z?i + hrn^ = 


-MF 


(S73) 


(S74) 
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This equation is a function of /i, which is a free parameter. Since it is a lower bound for the actual value of S, the 
best approximation would correspond to the value of h which maximizes the approximation. We have therefore that 
/i is a solution of the following equation 


f) f) f) f^Ttl 

0 = — log(SMF) = ^ log(^/t) + S'—(-/3m^ + hm) = S{h-2f3m) — 


(S75) 


where m is given by equation S73 We obtain therefore m = /i/(2/3) and then, by neglecting sub-leading terms in S 


and introducing m = h/(2j3) in equation S74 


- log Emf « log ( erfc(/i/2) exp 


By maximizing this equation respect to h we obtain 


d 


h f \ 


2 V/3 


d 


1 + /3 

4 p 


h fl 


0 = ^ log(^MF) = o ( a + 1 ) + ob log (erfc(/i/2)) = - ( - -h 1 ) - 


dh 




•y/7rerfc(/i/2) 


Equation S77 cannot be solved exactly. By expanding around h = 0 we obtain 


h fl 


which is solved by 


2 V/3 


h = 


1 


0=0 0+1 


2/3 


(S76) 


(S77) 


(S78) 


(S79) 


TT -I- /3(7r — 2) 

One can observe that the solution h = 0 corresponds to /3 = 0, i.e. to a non-interacting ecosystem. Expanding around 
h = 0 is therefore meaningful when the interactions are not too strong. It is possible to verify that the approximate 
solution |S79| is very close to the actual solution obtained by solving numerically equation |S77| also for not too small 
values of /3 


Using equation S79 into equation S76 we obtain 


— log ^mf 


P{l+P)Tr 


log erfc 


y+/3 


TT -I- P{tT — 2) 


(S80) 


(tt- f/3(7r - 2)) 

which is our final result. In figure [S^ we compare this equation with the volume computed numerically in the case of 
constant interactions, finding a very good match. 

In the most general case of an interaction matrix with nonconstant offdiagonal entries, we can consider equation |S72| 
as an approximation valid in the case of E 2 —>-0. As /3 was defined in terms of the generators, we can extend the 
approximation to the case 732 > 0 by considering /3 as the expected value of G’s entries, which corresponds to the 


average overlap of two rows of the interaction matrix {cos{r])), defined in equation SI16 In this more general case 
the mean-field value of S is expected to be a good approximation when var(cos(77)) is small enough. By substituting 
/3 = (cos(7y)), using equation S119 into equation S72| we obtain 

TrEi{2d - EiS) {2dEi + d-S {2El + E^)) 

(S81) 


S (d(2(7r - 2)£;i +f)-S {2{tt - l)El + TrE^)f 


log ( erfc 


y/^Ei{EiS - 2d) 


S (2(7r - l)Ef + TvEj) - d(2(7r - 2)£;i -f tt) 
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P 


Supplementary Figure S3: Approximation of H using mean field theory. The black dots are numerical simulations 


obtained by integrating H numerically (see section S41 for a constant interaction matrix. The red curve is the analytical 


approximation obtained using the mean-field approximation (see equation S811. /3 is a function of Ei and S, and is defined in 
equation 


S64 


The range of (3 considered here is the same of the one appearing in figure 1 of the main text. 


When var(cos(? 7 )) is not small, we observed that the empirical formula 


S 


log(^) 


TrEi{2d - EiS) {2dEi + d - S {2Ef + E^)) 
{d{2{TT - 2)Ei +7r)-S i2{TT-l)Ef+7rE^yf 

yE^Ei{EiS -2d) 


log erfc 


S (2(7r — i)Ef + ttE^) — d{2{TT — 2)Ei + tt) 


(S82) 


explains well the values obtained in simulations. This is the formula we used to make figure 2 in the main text. 

In order to simplify the expression and make it more readable, we can expand equation |S80| around /3 = 0, i.e., 


when the interactions between species are small. By expanding around /? = 0 and taking the logarithm of 

the expression, we obtain 


4log^MJ’ « log ( 1 - - 

b \ TT 


(S83) 


Equation 2 of the main text was obtained by substituting /? = (cos(7y)), using equation S119 in the case of E 2 = 0. 


S7. FEASIBILITY OF CONSUMER-RESOURCE COMMUNITIES 

This section considers explicitly a community with two trophic levels and consumer-resource interactions. While 
empirical communities have a more complicated interaction structure, this example is particularly relevant to better 
understand how S should be interpreted. 

We consider a system with Sr resource and Sc consumer (Sr + Sc = S) populations, whose dynamics is described 
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by equation SI with the interaction matrix 


f ■''1 ■ 

ZB^W 0 J 


(S84) 


where C is an Sn x Sn nonnegative matrix, B is an Sr x Sq nonnegative matrix, while Z and W are two positive 
diagonal matrices of dimension Sr x Sr and Sc Sc, respectively. 

If C is a positive diagonal matrix, any feasible hxed point is globally asymptotically stable [SS]. When C is not 
diagonal, one can prove that any feasible fixed point is globally asymptotically stable if CW~^ is positive definite 
(i.e., -CW-^ is negative dehnite). Assuming that this condition holds, stability of feasible fixed points is ensured 
and we can study feasibility alone. 

Using equation |S3[ we obtain the equations 

= E * + E ’ (S85) 


-rf = 

i=i 




(S86) 


where and are the intrinsic growth rates of resources and consumers, while n^* and are their equilibrium 
abundances. Since all the matrices that appear in this equation are nonnegative, an intrinsic growth rate vector is 
contained in the feasibility domain only if rf- > 0 for all i = 1,..., S'fl and rp < 0 for all i = 1,..., Sc- An intrinsic 
growth rate vector that does not respect these conditions is not in the feasibility domain. The feasibility domain is 
therefore fully contained in one orthant, implying that the maximum value of its size is S = 1. 

The 5'-dimensional volume of the feasibility domain is nonzero only if it is defined by S linearly independent 
generators. The generators of the feasibility domain are proportional to the columns of the interaction matrix. If 


the interaction matrix has the form of equation S84 Sr generators will have the form g = (t>,0), where v has Sc 
components. These generators can be linearly independent only if Sr > Sc, and therefore S > 0 only if Sr > Sc- 
More generally, if det(A) = 0, then S = 0 [57] . 

Assuming that the determinant of A is different from zero, we can use equation |S26| obtaining 
S = v/det(A) [ d^z expi^ZiAi 




(S87) 


'\i=l 


Given the structure of the matrix A, it is convenient to write z = (v,u), where v and u are two vectors with Sr and 
Sc components respectively. The argument of the exponential can be rewritten as 


E 


z^A.^Zj 


Sr Sr Sr Sc 

i — ~ E/ E/ ~ E/ E/ ~ 

i—1j—l i—1j—l 


(S88) 


By integrating over the variables it, we hnally obtain 


ViCijVj 




(S89) 
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Figure shows the size of the feasibility domain of a consumer-resource community, computed using Monte Carlo 
integration as explained in section |S4| We consider an interaction matrix with the structure of equation |S84[ with 
a diagonal C (i.e., C^- = 1 if z = j and zero otherwise) and scalar matrices Z and W (i.e., Za = Wu = rj). The 
elements of the rectangular matrix B were independently drawn from a lognormal distribution with mean fj, and 
variance c^/i^, where c„ is the coefficient of variation. Since C is equal to the identity matrix, then the interaction 
matrix is diagonally stable and therefore any feasible point is globally stable [SB]. Figure S4 shows the effect of r], /i 
and c-u on the size S of the feasibility domain. Interestingly, rj and ^ have a small effect on S, while the coefficient of 
variation has a strong influence on it. It is important to notice that, as explained above, as the interspecihc interaction 
goes to zero (and therefore both Cy and /r tend to zero), 2 —?► 0 as well. Note that in this case not all the species are 
self-regulated {An = 0 for the predators), and therefore in absence of interspecific interactions, 2 7^ 1. 


S8. EMPIRICAL NETWORKS AND RANDOMIZATIONS 

We considered 89 mutualistic networks and 15 food webs. Empirical networks are encoded in terms of adjacency 
matrices L, with = I if species j affects species i and zero otherwise. 


A. Mutualistic networks 


The 89 mutualistic networks (59 pollination networks and 30 seed-dispersal networks) were obtained from the 
Web of Life dataset (www.web-of-life.es), where references to the original works can be found. When the original 
network was not fully connected, we considered the largest connected component. 

In the case of mutualistic networks, the adjacency matrix L is bipartite, i.e., it has the structure 


L = 




(S90) 


where is a Sa x Sp matrix {Sa and Sp being the number of animals and plants respectively). The adjacency matrix 
contains information only about the interactions between animals and plants, but not about competition within plants 
or animals. 

We parameterized the interaction matrix in the following way: 


/ Lbo \ 

\ Ll o ) 


(S91) 


where the symbol o indicates the Hadamard or entrywise product (i.e., {AoB)ij = AijBij), while W^, 
and are all random matrices. and are both square matrices (of dimension Sa x Sa and Sp x S'p), 
while and are rectangular matrices of size Sax Sp and SpX Sa respectively. The diagonal elements 

and Wn were set to —I, while the pairs {W^, ^ji ) were drawn from a bivariate normal distribution 

with mean variance cr^ = cfiL, and correlation ptr^. Since these two matrices represent competitive interactions, 
fi- < 0. The the pairs {W ^^were extracted from a bivariate normal distribution with mean /r+, variance 
cr^ = cp,^, and correlation where fip > 0. 
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Supplementary Figure S4: Feasibility domain of consumer-resource community. We considered an interaction matrix of 
the form of equation |S84[ with Sr = 40 and Sc = 30, Cij = 1 if i = j and zero otherwise, ZiWj = r] > 0 for any i and j, and 
B with entries independently drawn from a Lognormal distribution with mean fi and variance For each parameterization 

we computed the feasibility domain H using the method explained in section [S^ The value of H is mostly determined by the 
coefficient of variation of the interaction, and it depends only weakly on the mean interaction strength jj, and the efhciency rj. 


We analyze more than 600 parameterizations, obtained by considering different values of /r_, /r+, c, and p. For 
each network and parametrization we computed the size of feasibility domain S. The bottom panel of Figure 2 in the 


main text was obtained by comparing obtained in this way with the analytical prediction obtained in equation S81 
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Supplementary Table SI: References and properties of the 15 food webs analyzed in the work 


Name 

S 

Number of links 

Connectance 

Ythan Estuary 

92 

414 

0.1 

St. Marks 

143 

1763 

0.17 

Grande Carigaie [60] 

163 

2048 

0.16 

Serengeti [61] 

170 

585 

0.04 

Flensburg Fjord [62] 

180 

1567 

0.1 

Otago Harbour |63| 

180 

1856 

0.12 

Little Rock Lake |64| 

181 

2316 

0.14 

Sylt tidal basin 

230 

3298 

0.12 

Caribbean Reef |66| 

249 

3293 

0.11 

Kongs Fjorden [57] 

270 

1632 

0.04 

Carpinteria Salt Marsh [55] 

273 

3878 

0.1 

San Quintin |68| 

290 

3934 

0.09 

Longh Hyne [55] 

349 

5088 

0.08 

Pnnta Banda [55] 

356 

5291 

0.09 

Weddell Sea [7D] 

488 

15435 

0.13 


B. Food webs 


A summary of the properties and reference of the food webs can be found in table |S1[ In the case of food webs the 
adjacency matrix L is not symmetric, and an entry Lij = 1 indicates that species j consumes species i. We removed 
all cannibalistic loops. Since both Lij and Lji are never simultaneously equal to one (there are no loops of length 
two), we parameterized the offdiagonal entries of A as 

= W+L,, + Wj-L,, , (S92) 

while the diagonal was fixed at —1. Both and W~ are random matrices, where the pairs {W^,W~j) are drawn 
from a bivariate normal distribution with marginal means and correlation matrix 

( i '’“f 'i (s») 

y cp± J 

We analyzed more than 200 parameterizations, obtained by considering different values of /x_, /i+, c, and p. For 
each network and parametrization we computed the size of feasibility domain S. The bottom panel of Figure 2 in the 
main text was obtained by comparing S obtained in this way with the analytical prediction obtained in equation|S81| 


In this case the analytical prediction overestimate the actual value of indicating that there is a role of structure in 
determining structural stability. 
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Supplementary Figure S5: Size of feasibility domain H in empirical networks and randomizations. Empirical networks and 
their randomizations were parametrized as explained in section [SS] Each empirical network was parametrized 100 times and 
the average H was compared with the one obtained by averaging 100 randomizations. Each point in this plot correspond 
therefore to a value of H of an empirical network and its randomizations averaged over the extraction of the interaction 
strenghts for a given combination of the parameters as explained in section 


S8 


S9. RANDOMIZATION OF EMPIRICAL NETWORKS: ASSESSING THE ROLE OF STRUCTURE 

A. Mutualistic networks 


We compared the size of the feasibility domain obtained for empirical networks with the corresponding random¬ 
izations. For each network we randomized the block 100 times, by generating connected networks with same size 
and number of links. We parameterized each randomized network independently as described in section [S8l and we 
compared their properties with those of the empirical network, parameterized independently 100 times. Figure |S5| 
shows the comparison between S of random and empirical networks. As expected from the fact that the analytical 
prediction for random matrices works well, the empirical values and the values obtained with randomizations are 
compatible. Comparing this figure with figure 2 of the main text we observe that the empirical values and the ones 
obtained with randomizations match also in the cases were the analytical approximation failed. This implies that 
the reason of the mismatch is due to the difference between the analytical approximation and the randomizations, 
and it is not due to the specific structure of the empirical interactions. There are two main sources of errors in this 
case. On one hand, ours analytical prediction is expected to work is the number of species is large enough and if the 
variance of interactions is not to high (that is not always true for the parametrizations used). On the other hand, 
our approximation was formulated for random matrices, while randomizations of mutualistic networks still conserve 
a bipartite structure. 

The randomization procedure explained above and figure |S5| show that the size of the coexistence domain obtained 
with empirical network structure is well predicted by the one obtained with random structure. This result does not 
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imply that structure has no effect on S, but it shows that, if this effect exists, it must be relatively small (compared 
for instance to the variation of S obtained by changing the interaction strengths), i.e. the relative error made by 
approximating empirical networks with random structure must be small. 

Since the effect of structure is small, it is also expected to be very sensible to the interaction strengths. When we 
parametrized empirical networks and their randomizations to obtain figure |S5[ we drawn the interaction strengths 
several times from a given distributions. The realized coefficients were therefore different across different networks, 
and the values of S shown in figure [S5| were averaged over these independent extractions. Since the difference between 
randomizations and empirical structure is small, it might be impossible to detect any difference with this procedure. 

In order to explore and quantify the effect of the empirical structure on the size of feasibility domain, we adopted a 
different parametrization and randomization method. Given an empirical network, we drawn the interaction strengths 


only once from a given distribution (as described in section S8). Using this list of interaction strenghts we parametrized 
100 times each empirical network. Different parametrization differ in the position of the coefficients, but not in their 
values that are conserved across parametrizations. We then compared their size of feasibility domain with the one 
obtained by parameterizing with the same list of coefficients 100 randomized networks obtained as explained above. 

Figures [S^ [S^ and [S^ show the results obtained for different distributions of interaction strengths (parametrized 
as explained in section [S^ . In absence of competition and in absence of variation in the interaction strengths, there is 
the maximum observable effect. As the competition level is increased and once variation in the interaction strengths 
is introduced, the effect of the network topology on the total size of feasibility domain becomes negligible. 


B. Food webs 

We compared the size of feasibility domain of empirical networks with their corresponding randomizations and a 
network generated accordingly to the cascade modeip?]. 

For each network, we randomized the adjacency matrix L 100 times, by generating connected networks with the 
same size and number of links. 

We also generated networks generated accordingly to the cascade model (using the same method explained in [Idjl. 
In this case the adjacency matrix was obtained by generating connected networks with the same size and number of 
links, by assigning a link between species i and j only if j > j. 

Figure [Sl0| is the same as figure 2 of the main text, with the addition of randomizations and networks generated 
with the cascade model. As expected the analytical prediction works very well in describing random networks, while it 
fails significantly to predict the size of the feasibility domain of cascade and empirical networks. To better quantity the 
difference between those empirical structures and randomizations, we compared each network separately in figure [STT] 
and |S12[ We observe that random networks have always larger feasibility domain than networks generated by the 
cascade model and the empirical ones. Networks generated via the cascade model almost always overestimate the 
empirical feasibility domains, showing that empirical network structure has a significant negative effect on the size of 
the feasibility domain. 
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Supplementary Figure S6: We measure the effect of mutualistic network structure on the size of the feasibility domain as 
described in section [S9A[ Red violin plots are randomizations, green ones are empirical networks. The empirical networks are 
grouped in four rows based on the number of species {S < 50, 50 < S < 80, 80 < 5 < 150 and S > 150, respectively). This 
figure was obtained with = 0.25fj,max, A*- — 0 s-nd for three different values of c. 
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Supplementary Figure S7: Same as figure 


S6 


but with jj.+ = 0.25jj,max and /i_ = 0.5/i+ 
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Supplementary Figure S8: Same as figure 


S6 


but with — 0.5fj.max and fj— = 0 
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Supplementary Figure S9: Same as figure 


S6 


but with = 0.5fimax = M- = 0.5^+ 
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# Empirical 

# Randomized 

# Cascade 


Supplementary Figure SIO: In this figure we compared the analytical prediction of the feasibility domain obtained in 
section [S6| with the numerical calculated values for random networks, empirical networks and networks generated via the 
cascade models. The feasibility domain of random networks is well predicted by our analytical approximation, which fails to 
predict the empirical one and the one obtained using the cascade model. 

SIO. POSSIBLE BIASES IN PREVIOUS ANALYSIS OF STRUCTURAL STABILITY 


In section |S4| we showed how to estimate the feasibility domain numerically in a fast and reliable way. In previous 
approaches [S], the feasibility domain (structural stability) was not directly calculated, but approximately inferred 
using a regression method. In this section we show that the method used by Rohr et al. [9] could be biased and is 
not always applicable. 

The Authors considered a bipartite mutualistic system described by the dynamical model 


dnf 

dt 

dnf 

dt 


= rii 


Sa 

i=i 

Sp 


Punf 


p p 






nj '^J 

'yP-n^ 

’^J J 


(S94) 


1 + '■f it’d 

where Sa (Sp) is the number of animals (plants), and nf (nP) is the abundance of animal (plant) species i. For the 
purposes of this section we consider the case of linear functional responses hf = hf = 0, as all the methodology used 






46 


0.0 

- 0.2 

-0.4 

- 0.6 

- 0.8 

0.0 

- 0.2 

-0.4 

- 0.6 

0.1 

0.0 

- 0.1 

- 0.2 

-0.3 


[■) 


D) 

O 












4 


4 


1 




-4 




0 


❖ 




❖ 




6o 


4 


6 


6 


Empirical Randomized [^Cascade 


Supplementary Figure Sll: We measure the effect of food web network structure on the size of the feasibility domain as 
described in section [S9 B[ Red violin plots are randomizations, green ones are empirical networks, while blue ones correspond 
to the cascade model. This figure was obtained as explained in section fSQ B| with /i_ = 0.25^max, IJ,+ = 0.5fi- and for three 

different values of c. 
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in [S] was developed in this case. If the functional response 
interaction matrix A is given by 


is linear, this equation reduces to equation SI 


where the 


A = 




(S95) 


Here 13^ and /3^ are Sa x Sa and Sp x Sp matrices, respectively, while 7"^ and 7^ are Sa x Sp and Sp x Sa 
matrices. The Authors used a constant parameterization for the competition parameters, setting /id = //T = l and 
/3d = = p li j ^ i. The mutualistic benefits were parameterized as 


.-yd 


= 7o 


Tin = 70 


{fY 

{kfY 


(S96) 


where Ljy is the nonzero block of the adjacency matrix of the interaction network, i.e., Ly- = 1 if there is an interaction 
between animal i and plant j, and zero otherwise. The numbers kf = ^ij kf = Ylj=i degree 

of animal/plant i. The two remaining parameters, 70 and <5, quantify the levels of mutualistic strength and the 
mutualistic tradeoff m- 

The method proposed by Rohr et al. [S] was based on what the Authors called the “structural vector”. It was defined 
as the center of feasibility domain and was calculated by transforming the mutualistic dynamics into an effective 
competitive one. Using this effective dynamics it was possible to calculate an effective structural vector, which was 
then transformed back to the one of the mutualistic system. Starting from the structural vector, the Authors considered 
different perturbations of the growth rates by changing their direction from that of the original structural vector by 
some given angle. The dynamics was then integrated and the probability that all species survived was calculated, given 
a particular perturbation. Running this across several different perturbations and parameterizations, it was possible 
to perform a regression between the interaction parameters, the angle by which the growth rates were perturbed, 
nestedness, and other parameters appearing in the interaction matrix. Using the coefficients obtained through the 
regression, it was quantified the effect of nestedness and other properties on the size of the feasibility domain. Here 
we present some possible issues emerging from this approach. 

It is not al-ways possible to find the structural vector. In order to calculate the structural vector, one needs 
to transform the mutualistic system into an effective competitive one. One can define the matrix T = 1 + 7^”^, 
where 1 is the identity matrix and 

(3 = 

and 

7 = 

By multiplying both sides of equation |S95| by T one obtains the effective interaction matrix 


/3^ 0 
0 /3^ 


7 


0 7"^ 

P 0 


(S97) 


(S98) 


Acff — 


-/3^+7'^(/3'^)-V 0 

0 -/3^+7^(/3^)-i7^ 


1 

to 

0 1 

1 

0 

sfff / 


(S99) 
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In order to calculate the structural vectors, one has to assume that the eigenvectors associated with the largest 
singular eigenvalues of and have only positive components (and an equivalent condition on 

Bgg ). This is not generally true, as also stated by the Authors [S]. They therefore imposed the extra assumption that 
and B^g{B^g)'^ indeed have only positive entries (and the equivalent conditions on B^g). In this case, 
the Perron-Frobenius theorem allows all entries of the leading eigenvector to be chosen positive; i.e., it necessarily 
points in some feasible direction. The Authors then identified the structural vectors with these eigenvectors. 

However, the extra requirement that {B^g)"^B^g and B^g{B^g)'^ be strictly positive imposes constraints on the 
interaction matrix that reduces the number of parameterizations that can be analyzed with this method. Since this 
assumption does not hold in general, there are cases in which the structural vector does not exist. Using our approach. 


this vector is not needed (see sections S4 and Sll). 

When the structural vector exists, it is not unique. Under what conditions would the matrices 
and satisfy the conditions of the Perron-Frobenius theorem? It is easy to show that this can never be the 

we see that A®® is block-diagonal, therefore (A®®)^A®® and A®®(A®®)^ are block-diagonal 


S99 


case. From equation 

as well. This means that the Perron-Frobenius theorem does not hold (the matrix is reducible); instead, the two 
diagonal blocks each have an all-positive leading eigenvector (assuming that all the coefficients are positive in the two 
blocks). Any linear combination of the two will have positive components. There is no reason to prefer one linear 
combination over another, and while it is true that some linear combinations may point closer to the center of the 
feasibility domain, there is no way to determine using the Authors’ methods which combination does, if any. 

The structural vector is uot the ceuter of the feasibility doruaiu. Let us assume now that the structural 
vector exists and it points toward the center of the feasibility domain of the effective competitive system. To obtain 
the structural vector, one has to transform it back to a vector of the original, mutualistic system. The transformation 
from the effective to the original system is done by multiplying with the matrix This matrix is not a rotation, 

and therefore it does not preserve the angles between vectors. Even if a vector is the center of the feasibility domain in 
the effective system, it will not in general be the center of the original domain. In particular, its distance to the actual 
center of the original domain will be dependent on parameterization and network structure, as the transformation 
matrix depends on these. 

In contrast, the center of the feasibility domain can be easily expressed with our approach in terms of the matrix 
A and its associated generators (section S3 equation |S20[). It is also easy to check that the barycenter is different 


from the one obtained using the method of Rohr et al. [5]. 

The regressiou procedure cau iu priuciple produce biases. The relationship between network structure 
and the size of the feasibility domain was obtained by calculating the probability of coexistence p( 0 a,^p), where 
(^A/p is the angle by which the direction of the growth rate vector of animals/plants was changed with respect to the 
structural vector. The Authors then performed a linear regression 


logit(p( 6 »A, dp)) ^ /3i log 9a + I32 log 9b + kloC + 

+ PbloCN + + PsjCS'^ 


(SlOO) 


where C is the connectance of the mutualistic adjacency matrix and N is its nestedness (note that 7 , used by Rohr 
et al. [5], is equal to Cyo). The fitted parameters where then used to determine the effect of nestedness and other 
quantities on the feasibility domain. The functional dependence assumed above cannot be justified a priori, and an 
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incorrect functional dependence can in principle lead to erroneous fitting results. For instance, the effect of those 
properties could be different depending not just on the raw angle of perturbation, but also which direction that angle 
is taken in. We can imagine two feasibility regions with the exact same size but different shapes: one of the two is 
equally wide in all directions, while the other stretches very wide in some directions but is extremely narrow in others 


(see section Sll). For sufficiently small values of dyi/p, one will never leave the feasible domain in the first of these 
examples, but may do so in the second if the perturbation is performed in one of the “narrow” directions. The first of 
these cases will therefore appear more feasible than the second, even though the total size of the two feasibility regions 
is in fact the same. On the other hand, if the values of 6aip s^re large enough, than the perturbed vector in the first 
case will never be feasible, while it will be feasible in the second case because of the “wide” directions. Moreover, this 
method does not allow one to calculate the feasibility domain for a given network and parameterization, as one can 
calculate only the probability of coexistence given an angle of perturbation. 


Sll. DISTRIBUTION OF SIDE LENGTHS 


In section |S3| we showed that the feasibility domain is a convex polyhedral cone in the space of intrinsic growth 
rates r. Since the stationary solution of equation |S1| is linear in r, we can study the feasibility domain considering 


only vectors on the unit sphere’s surface. In section S4 we dehned S, which quantifies the volume of the feasibility 
domain. 

The size of the feasibility domain, i.e., how many combinations of the intrinsic growth rates correspond to a feasible 
fixed point, is not the only interesting property. Two systems having the same number of feasible combinations 
of growth rates (i.e., the same value of S), can respond very differently to perturbations of the growth rates. We 
imagine here that a perturbation (e.g., a change of the abiotic conditions) correspond to a change in the growth rate 
vector. Since we can consider normalized growth rate vectors (because of the linearity of the equations), the effect of 
a perturbation on feasibility depends only on the angular change of the growth rate vector and not on its length. 

The volume S quantifies how many growth rate vectors are compatible with coexistence. Let us consider a feasible 
growth rate vector, and perturb it in a random direction. What is the probability that the new vector is still feasible? 
This is not just a function of the size S of the feasibility domain. Indeed, one can imagine that the feasibility domain 
is about equally spread in every direction—or that, for the exact same value of S, the feasibility domain is streched 
in some directions but is very narrow in some other ones. A perturbation in one of the “narrow” directions is much 
more likely to lead out of the feasibility domain in the latter case than in the former. 

To quantify this property, one strategy could be to measure the different responses on the perturbation (i.e., the 
probability of being feasible) depending on the direction of the perturbation (in which direction we change the growth 
rate vector). This choice has the big disadvantage of depending not only on the properties of interactions (the 
interaction matrix A ), but also on the strength of the perturbation (the angular displacement between the initial and 
the final growth rate vector) and the growth rate vector before the perturbation (e.g., if the initial vector is close or 
far from the edge of the feasibility domain). We propose instead a purely geometrical method to quantify the response 
to different perturbations (see figure 1 of the main text). 

The feasibility domain, when restricted to the surface of a hypersphere, can be imagined as the generalization of a 


triangle on a sphere (see section S13). The natural, geometric quantities bounding the maximal perturbation that will 
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leave the system feasible, are the lengths of the triangle’s sides. When S species are considered, there are S{S — l)/2 
sides. Their lengths measure the maximum permissible perturbation of the growth rates in the corresponding direction 
if one is to retain feasibility. This property has the advantage of being purely geometrical, depending only on the 
interactions (via the interaction matrix) and not, for instance, on any choice of the initial conditions. 

We can measure the distribution of the side lengths. Imagine we have two interaction matrices with the same 
S, but with very different distributions of side lengths. One of them has all sides of equal length, while the other 
one has a more heterogeneous distribution. In the first case any direction of the perturbation is expected to have a 
similar effect, and there are no particularly dangerous directions. In the second case there are some directions of the 
perturbation that are much more dangerous than others, and even a small change of conditions along one of those 
dangerous direction can lead to the extinction of one or more species. 


We know that the feasibility domain is a convex polyhedral cone (see section S3). Its “corners” are identified by 


its generators and its sides are determined by all pairs of generators (see section S13 for the S' = 3 case). 

Since we are considering growth rates on the unit (hyper)sphere, and the generators are normalized to one, any pair 
of generators will lie on the sphere’s surface. The scalar product of two generators is the cosine of the angle between 
the two. Since the two generators are on the unit ball’s surface, the arc between the two (which is the side length) is 
equal to the angle. We have therefore that the length of the side of the feasibility domain corresponding to a pair of 
generators g® and is 


gij = arccos (g* • g-’) . 


(SlOl) 


Using equation S12 we can express the S{S — l)/2 side lengths of the convex polytope explicitly in terms of the 


interaction matrix: 


rjij = arccos 


^kiA-ki 


(S102) 


We are interested in the distribution of the side lengths, and in particular in its heterogeneity. In the following section 
we will calculate these quantities for random matrices. 


A. The distribution of side lengths in random matrices 


In this section we obtain the distribution of sides length for large random matrices, whose entries are distributed 
accordingly to an arbitrarily bivariate distribution. 

We assume that the diagonal elements of A are all equal to —d (this hypothesis can be easily generalized), while 
the offdiagonal pairs (Aij,Aji) are random variables with distribution q{x,y). Our goal is to find the distribution of 
the side lengths g in the large S limit, defined as 


P{il) = c; rdrE/ n( dA 

mn dA 

nm^(-^mn7 -^n 

m>n 

-^ki-^kj \ \ 


S^oo S{S-l) 
X d g — arccos 


^kiAki ^ 


(S103) 


Since we are summing over all i and j, and all the rows are identically distributed, we can remove the sum and 
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consider just two rows: 


P{v) = Jus dA 

mn d^ 

nm qiA 

mm 

J2k ^klAk2 


m^n 


(S104) 


X (5 ?7 — arccos 


_ \/J2k ^klAkl J2l ^12^12 ^ 
Since we are interested in the large S limit, we have that 


AkiAki =Aii + Ki -d+ {S - 1) / dxdy q{x, y) 

k k>i 


(S105) 


= -d+{S-l){El+El), 

where Ei and E 2 are the first and second marginal moments of q (equations |S35 and S36). Let us call this quantity 


Z. In this limit we therefore obtain 


J2k ■^klAk2 


P{jl) —^ iv arccos 

m>n 

= j^i^ / n [<^^mndAnmqiAmn, Anmij Z |sin(?7)| 6 I Zc0s(77) - AklAk2 

m>n \ k 

= Z\ sin(77)| ^1^ / (^dAmndAnmq{Amn, Anm)j S I Z cos(ry) - ^AklAk2 

m>n \ k / 

= ,Z^| siii(7y) I 

m>n 

X 6 iz COs{r]) — ^11^21 — ^22^12 — AklAk2 I 

V k>2 J 

= ^|sin(/y)| (^dAYnndAnmq{Amn-) Anm)^ 

X S izcos{r]) + d{Ai2 + A21) — AkiAk2 j 

= Z|sin(?7)| J dt J ds J dAi2dA2iq{Ai2, A2i)S{t - A12 - A21) 

X / 77 dAkidAk2qiAki)q{Ak2)S I s — AkiAk2 j 

k>2 V fc>2 / 

X 6 iz cos{r]) + dt — E Xlfel^fe2 I 

V k>2 J 

= Z|sin(?7)| J dt J ds J dxdy q{x,y)S{t - {x + y)) 

X y ( 77 dzkdwkq{zk)q{wk) I ^ I s " E ZkWk j d{Z cos(?7) + dt — s) , 


vfc=l / 

where q{z) is the marginal distribution of q{x,y): 


,c)=/dx,(x,.)=y'dx 


We can introduce the distribution of the sum: 


qs{t) = J dxd yq{x, y)S{t - {x + y)). 


(S106) 


(S107) 


(S108) 
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The term 


/S-2 \ / S-2 

dzfcduifc q{zk)q{wk) <5 s - 

^fc=l / V k=l 


ZkWk 


(S109) 


is the distribution of a sum of S' — 2 uncorrelated random variables. These random variables are the product zw of 
two random variables whose distribution is q. Since the second moment of q{x) is finite, the central limit theorem 
holds and this distribution converges, in the large S limit, to a Gaussian distribution with mean 


and variance 


We have therefore 


P{r]) = Z|sin(77)| J dtds qs{t) 


S J dxdy q{y)q{x) xy = SEf 
S (^J dxdy q{y)q{x) {xyf - EI^ = SE^ . 

V ■^SE^ ) 


(SllO) 


(Sill) 




5{Z cos(? 7 ) + dt — s) = {S(EI + E 2 ) — d) 


|sin(? 7 )| 


f ( {S{El+El)cos{ri)-dcos{ri)-SEl + dty 

j “P-2Sl- 


(S112) 


The distribution of 77 is not universal as it depends on qs(t), which depends on the distribution of the coefficients. On 
the other hand, the dependence is explicit, and it is possible to calculate P{r]) for any distribution q{x,y). 

We show explicitly the case of q(x,y) being a bivariate normal distribution, i.e., 


q[x,y) = 


1 




exp - 


{x - Elf + {y- Elf - 2Efx - Eify - Ef 


2El 


(S113) 


In this case qsf) is a normal distribution, and can be obtained from eq S108 

{t-y- Elf P{y- Elf - 2Eft - y - Ei){y - Ef 


il-Ef{t-2Eif 
4Ei 


2El 


= exp - 


1 


(S114) 


2v^S2(l + S,)ynr^ ■ 


Substituting into equation S112 we see that P{rj) has the form of a convolution of two Gaussians, and turns out to 
be equal to 

'cos{ri) - (cos(? 7 )))^\ 


Pfl) = 


■\/^var(cos(? 7 )) 


exp 


2 var(cos(77)) 


(S115) 


The mean {cosirf}) and variance var(cos(? 7 )) will be computed in the next section in the most general case of an 
arbitrary interaction distribution. 


B. Moments for random matrices 

As explained in the previous section, the distribution of the side lengths is not a universal quantity, as it depends 
on the distribution of interaction strengths. In this section we compute the mean and the variance in the general case, 
showing that they depends only on Ei, E 2 and Ec- 
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Here and in the main text we do not report the moments of the side length 77 , but the moments of its cosine. The 
cosine of the side length measures the overlap between two rows of the interaction matrix (or the scalar product of 
two generators of the convex polytope). As its value gets close to one, the side length approaches zero. 

Starting from equation S102[ we have that 

-^ikAjk \ 


(cos( 77 )) = 


1 

S{S-1) 


^cos(77*j) = 






Since we are interested in the large S limit, we can write the denominator as in equation |S105| and obtain 

1 

S{S-l) 




i¥=j 


-d + SiEf+El)J ’ 


and then 


(cos( 77 )) = 


In the large S limit, this becomes 


1 


S{S-1) 


E 


^uAji + AijAjj + ^ikAj 


H'fc 


{cos{r])) = 


-d + S{El + El) 


-2dEi + SEl 


-d+iS-2){Ef + El) 

to leading order in S. 

In a similar way, we can write the second moment as 


( cos( 77)2) = ^cos(r7,,)2 = ^ 


^ikAjk 




1 ) \ V^k ^ikAik , 


In the large S limit we obtain 
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We can compute the averages of the different terms, obtaining 

1 . . .0 1 
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+ Aji)"^ — E^^b + + 2AjiAj,) 
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(S116) 


(S117) 


(S118) 


(S119) 


(S120) 


(S121) 


(S122) 
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(S123) 
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and 


5'(5'_ 1) X! ( X! 


jk 


We finally get that, in the large S limit, 
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X! X! X! AikAuAjkA^ 
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(S124) 


(S125) 


S12. SIDE HETEROGENEITY FOR DIFFERENT STRUCTURES AND EMPIRICAL NETWORKS 


In figure [ST^ we considered the effect of four nonrandom structures on the mean and variance of the side lengths. 
The interaction strengths were drawn from a normal distribution with given mean, variance, and correlation. For some 
structures we considered multiple interaction types and therefore multiple means (one positive and one negative), in 
which case the coefficient of variation of the interactions and the correlation was constant and independent of the 
mean. Networks were parametrized as explained in section [S8l 

• Modular. In this case we considered interaction matrices with a perfect block structure (to generate figure 3 
we considered four blocks of equal size). 

• Bipartite. In this case we considered an interaction matrix with two bipartite blocks of equal size. The mean 
interaction of the offdiagonal blocks was set to be negative, while the one of the in-diagonal blocks was positive. 

• Nested. The interaction matrix had a bipartite structure. The diagonal blocks had a random structure with 
negative mean interaction strength. In the offdiagonal blocks, we consider a connectance equal to one half and 
we built a perfectly nested matrix. The mean interaction strength was positive in the offdiagonal blocks. 

• Cascade. We build a matrix using the cascade model, and parameterize it with a positive and a negative mean 
depending on the role of the species in the interaction. 


In the case of empirical structures, figure 3 of the main text, was obtained considering the same networks and the 


same parameterizations considered in section S8 We compared var(cos(77)) with the values expected in the random 
case. Figure 


S14 


shows the comparison between ( cos{ri f obtained for empirical networks with the null prediction. 
Its value is well predicted by the null expectation for mutualistic networks, while the null expectations underestimates 
this value for food webs. This is consistent with the fact that the size of feasibility domain of random networks is 
larger that the one of empirical networks. 


S13. FEASIBILITY DOMAIN FOR S = 3 

When S' = 3, it is possible to visualize in three dimensions a convex polyhedral cone and the feasibility domain [25) . 
In figure [STS] we show a convex polyhedral cone in three dimensions and its generators. 
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Null <cos(ri)> 



0.00 0.02 0.04 0.06 0.08 


0.00 0.02 0.04 0.06 0.08 



0.000 0.025 0.050 0.075 0.100 0.00 0.02 0.04 0.06 0.08 

Null fluctuations of cos(ri) 


Supplementary Figure S13: We measure the effect of non random structures of mean and variance of side lengths. With 
the exception of the cascade model, all the structures considered do not have an important effect on (cos(?7)). On the other 
side, all the non-random structures considered have a positive effect on the fluctuations of cos{rj). All the networks considered 

had a connectance C = 0 . 2 . 


An important feature of convex polyhedral cones is that if r belongs to the cone, then so does cr for any positive 


constant c. As explained in section [S3l this is a consequence of the linearity of equation SI It is relevant therefore to 
limit our analysis to the growth rate vectors on the unit sphere, i.e., to vectors r such that 


r = 


= 1 . 


(S126) 


When we consider the vector in the feasibility domain on the surface of a unit sphere we obtain the areas of figure 1 
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Food webs Mutualistic 




Supplementary Figure S14: Comparision between (cos(» 7 )) obtained for empirical networks and its null expectation for 
empirical food webs and mutualistic networks. This figure was realized with the same parametrization of figure 3 of the main 

text and as described in section [i 


S8 



Supplementary Figure S15: Convex polyhedral cone and its section on a sphere. Left: the feasibility domain is a convex 
polyhedral cone, which is completely determined by its S generators (when 5 = 3 we have 3 generators , and g^). 

Center: since we consider a linear equation we can focus the analysis only on the intersection between the convex polyhedral 
cone and the unit sphere’s surface, which in three dimensions results in a spherical triangle. Right: each side of the convex 
polyhedral cone can be determined from a pair of generators as an arc rj of the sphere’s surface. Since we are considering the 
unit sphere, the arc length g is equal to the angle between the two generators. 
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in the main text. In this case, the quantity S is the area of the triangle, while the side lengths are the three sides of 
the triangle. Note that the polygon is not a triangle (as it lies on a sphere), but rather a spherical triangle. Its sides 
are arcs of a circumference, while its corners are identified by the three generators of the convex polyhedral cone. 

In the S = 3 case it is possible to obtain a closed expression for the area S [46] : 


„ 8 

.n = — arctan 

TT 


( 


|det(G)| 

^ + 





(S127) 


where the second term adds one to the first term when the argument of the arctangent is negative, while the matrix 
G is defined as 


Gu = g^ ■ g^ . 


Equation SI27 can be expressed directly in terms of the matrix A using equation SI2 


(S128) 


S14. NONLINEAR PER CAPITA GROWTH RATES 


In general, the effect of a species on the per capita growth rate of other species is not linear. Equation [ST] assumes 
this to be linear and the results presented in this paper were obtained under this assumption. Nonlinearity of the per 
capita growth rates can be thought of as a dependence of the interaction matrix A on n: 


drq 

dt 


i=i 


(S129) 


For instance, in the case of predator-prey interactions with a Holling type II functional response, it would have the 
form 


Aij{n) = 




1 + 


(S130) 


where the hij are the handling times. 

The presence of nonlinearity has strong consequences for both feasibility and stability. It is no longer possible to 
disentangle feasibility and stability with a simple condition on This means that feasibility will depend not only 
on the direction of r, but also on its length. 

The results presented here are a necessary stepping stone for assessing the feasibility of nonlinear systems. When 
the degree of nonlinearity is small (e.g., hij « 0), one can use our results, valid for the case hij = 0, to find the center 
of the feasibility domain and the generators. One can then treat the departure from hij = 0 as a small perturbation, 
and therefore, instead of having to explore the full vast parameter space, use the solution of the linear case as a 
starting point for numerical calculations to converge on the actual, nonlinear feasibility domain. On the other hand. 


in the limit of very large hij values. It is possible to show that the nonlinear form in equation S130 is approximately 


linear, and so again it is possible to use our method. The effect of intermediate values of hij on the feasibility domain 
is, however, still an open question. 
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